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ABSTRACT 

This  report  gives  detailed  descriptions  and  instructions  for  use  of 
a  system  of  programs  for  crystallographic  calculations,  including  least- 
squares  refinement  with  generalized  systems  of  constraints,  calculation 
of  bond  distances  and  angles  with  errors,  Fourier  synthesis,  plotting 
of  contours  in  Fourier  maps,  and  preparation  of  structure  factor  tables 
for  publication. 
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These  programs,  although  tested  on  many  problems  and  several 
computing  systems,  may  not  give  correct  results  in  all  circumstances. 
Therefore,  neither  the  authors  nor  their  institutions  will  guarantee  the 
correctness  of  the  program  in  all  respects,  but  the  authors  will  greatly 
appreciate  being  informed  of  discrepancies  in  the  code  or  desirable 
'options  not  currently  included. 
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RFINE4:  A  program  for  least  squares  refinement 
including  general  constraints 


General  Description 

This  program,  called  RFINE4,  performs  the  calculations  necessary 
to  refine  crystal  structures  using  the  techniques  of  full-matrix 
least-squares  refinement.  In  addition,  it  can  compute  bond  distances 
and  angles  or  thermal  ellipsoids  so  that  these  functions  may  be 
monitored  throughout  the  various  stages  of  a  refinement.   File 
control  capabilities  have  been  designed  s'o  that  parameters  may 
be  saved  at  the  end  of  a  given  computer  run  or,  alternatively,  virtually 
any  parameter  may  be  changed  at  almost  any  time. 

A  scheme  of  derivative  modification  and  parameter  restoration 
which  eliminates  the  need  for  coding  special  subroutines  to  treat 
constrained  parameters  has  been  included  in  this  program.  This  method 
will  handle  thermal  and  positional  restrictions  for  atoms  in  positions 
of  special  symmetry  and,  in  conjunction  with  the  occupancy  refinement 
capability,  a  multiposition,  partially  ordered  structure  may  be  con- 
strained to  a  particular  chemical  composition.  Provision  has  been 
made  for  incorporating  parameters  which  are  not  members  of  the  set 
used  in  conventional  least-squares  refinement .  Any  set  of  linear 
constraints  may  be  specified,  and  also  groups  of  atoms  may  be  treated 
as  rigid  bodies  both  with  respect  to  position  and  orientation  (shape 
constraints)  and  with  respect  to  thermal  motion. 

In  its  present  form,  the  program  allows  a  maximum  of  60  atoms, 
20  chemical  species ,  10  scale  factors ,  1320  parameter  dependencies 
and  200  refined  parameters  for  a  sfructure  in  any  space  group.  With 
the  maximum  number  of  parameters ' and  no  use  of  overlays  the  program 
requires  approximately  65000  words  of  core  storage.  However  more  restric- 
ted problems  can  be  operated  in  smaller  core  sizes,  and  the  program 
can  be  further  reduced  in  size  by  means  of  overlays.  A  suggested  overlay 
structure  is  given  in  Appendix  A.   Some  operating  systems  allow  for 
dynamic  allocation  of  core.  Such  provisions  can  be  used  to  advantage. 


File  and  Record  Structure 

In  the  following  discussion  a  record  is  the  unit  of  information 
formally  defined  as  a  logical  record,  i.e.  the  amount  of  information 
transferred  in  a  single  non-formatted  read  or  write  statement.  This 
description  is  independent  of  any  physical  device  which  may  be  associated 
with  the  file  and  the  program  has  been  coded  so  as  to  operate  fairly 
efficiently  with  any  type  of  auxiliary  storage  currently  available. 
Furthermore,  a  data  set  shall  mean  the  collection  of  records  needed 
to  describe  a  structure.  Any  reference  to  file  or  data  labels  should 
not  be  confused  with  file  control  information  for  a  particular  computer 
system. 


Master  File.   The  file  referred  to  as  ISC1  and  assigned  to  logical 
unit  10  is  the  main  file  of  the  program.   This  file  consists  of  a  data 
set  containing  all  the  information  required  to  resume  least-squares 
refinement .  This  data  set  has  two  label  records ,  one  or  more  structure 
factor  records,  a  constraint  relations  record,  a  variance-covariance 
record,  and  a  record  containing  variances  of  previously  refined  parameters 
(useful  in  large  problems  where  the  refinement  must  be  "blocked") . 
The  label  records  contain  a  12-character  alphanumeric  identification 
tag,  a  cycle  number  and  all  the  atom  parameters,  symmetry  data  and 
scattering  factor  information  needed  for  this  structure.  The  reflection 
data  consist  of  the  Miller  indices,  the  observed  structure  factor  on 
a  relative  scale,  the  standard  deviation  of  F0bs»  tne  value  of  sin28/A2 
for  the  reflection,  an  integer  specifying  which  scale  factor  applies 
to  this  reflection,  3 — the  extinction  parameter,  a  flag  to  mark  reflections 
which  had  an  intensity  less  than  the  minimum  observable,  the  dif fractometer 
angles  x  ,  $  and  u  +  6  in  radians ,  and  the  value  of  each  of  the  scattering 
curves  for  this  reflection.  These  quantities  are  packed  into  a  buffer, 
starting  with  word  2,  until  another  reflection  would  require  the  buffer 
to  be  larger  than  1023  words.  An  end  of  reflection  records  marker, 
which  is  non-zero  only  in  the  last  block,  is  contained  in  word  1. 
The  information  for  a  single  reflection  is  never  split  between  two 
records  and  all  records,  including  the  last  which  may  be  only  partially 
filled,  are  the  same  length.  The  end  of  the  reflection  list  is  marked 
by  a  Miller  index  h  having  a  value  of  3000.  If  the  number  of  scattering 
curves  is  N,  the  number  of  reflections  per  record  (n)  may  be  found 
from  the  integral  part  of  1022/ (N  +  12).  The  length  of  the  record 
is  then  1023  minus  1022  modulo  n.  The  variance-covariance  record 
contains  the  variance-covariance  matrix  for  the  previous  cycle  of  refinement 
in  packed,  upper  triangular  format. 

Intermediate  File.  This  file,  named  ISC2  and  assigned  to  logical 
unit  29,  contains  the  structure  factor  records  for  the  current  structure. 
A  distinct  file  allows  the  scattering  factor  information  to  be  changed 
without  reading  the  structure  factors  from  cards. 

Auxiliary  Files.   A  file  named  ISC3  and  assigned  to  logical  unit 
11  is  used  as  an  alternative  structure  factor  input.  A  file  named 
ISC4  and  assigned  to  logical  unit  12  is  used  to  pass  information  to 
Fouier  synthesis,  structure  factor  tabulation  and  statistical  analysis 
programs.  A  file  named  ISC5  and  assigned  to  logical  unit  28  is  used 
for  temporary  storage  of  constraint  relationships. 

System  Files.   This  program  uses  a  file  named  IN,  logical  unit 
5,  as  a  card  read  file,  a  file  named  IOUT,  logical  unit  6,  as  the  print 
file,  and  a  file  named  IPUN,  logical  unit  1,  as  the  punch  file.  These 
assignments  may  be  changed  for  a  particular  system  by  a  block  data 
subprogram. 

Description  of  Routines 

Main  Program.   The  main  program  serves  as  a  driver  for  the  least- 
squares  section.   It  reads  one  card  which  contains  a  12-character  label. 


If  a  new  data  set  is  being  constructed,  this  label  will  be  inserted 
into  that  file.  Control  is  then  transferred  to  subroutine  RFINE.   If 
the  updated  information  is  written  back  to  file  ISC1,  the  subroutine 
exits  normally.   In  that  event,  the  main  program  prints  a  message  stating 
that  the  file  has  been  updated. 

Subroutine  RFINE .   This  routine  is  the  primary  one  for  the  structure 
factor^and  least-squares  calculation  and  is  responsible  for  controlling 
the  functions  computed  and  the  parameters  used . 

When  this  routine  is  first  entered,  initialization  is  accomplished 
and  the  main  cycle  loop  is  entered.  The  first  step  is  to  read  the  title 
and  control  cards  for  the  cycle.   If  this  is  the  first  cycle  of  the 
current  run,  the  new  parameter  indicators  from  the  control  card  are 
interrogated  and  if  all  are  non-zero,  signifying  that  all  parameters 
are  to  be  loaded  from  cards ,  the  data  set  on  file  ISC1  is  presumed  not 
to  exist;  however,  if  any  of  the  indicators  are  zero,  the  first  two  records 
of  that  data  set  are  loaded  into  memory.  Next,  any  new  parameters  required 
are  read  from  cards.  As  new  parameters  are  read,  the  values  are  printed. 

If  the  unit  cell  parameters  are  to  be  changed,  that  card  is  read 
and  the  direct  space  metric  tensor  is  prepared.  This  is  a  3  x  3  symmetric 
matrix  whose  elements  are  the  scalar  products  of  the  unit  cell  vectors. 
This  matrix  is  then  inverted  using  subroutine  SYMINV  described  below. 
The  inverse  matrix  is  then  the  reciprocal  space  metric  tensor  and  the 
determinant  of  the  original  matrix  is  the  square  of  the  unit  cell  volume. 
Note:  If  the  unit  cell  parameters  are  to  be  changed,  the  structure  factors 
must  also  be  changed  if  the  effects  of  the  cell  data  are  to  be  applied 
to  the  values  of  sin  0/X  . 

New  symmetry  cards,  if  required,  are  now  read.   Each  card  lists 
an  equivalent  position  in  very  nearly  the  same  form  as  Vol.  I  of  the 
International  Tables  with  only  one  major  difference;  positions  related 
by  a  center  of  inversion  or  lattice  centering  must  not  be  included  since 
they  are  internally  generated.  Using  this  scheme,  all  space  groups 
may  be  generated  by  using  a  maximum  of  24  cards.   Each  symmetry  data 
card  is  transformed  into  matrices  R  and  T  satisfying  the  equation 

X'  =  RX  +  T 

with  X1  equal  to  the  vector  of  transformed  coordinates,  X  is  the  position 

vector  with  components,  xyz,  R  is  a  3  x  3  rotational  matrix  and  T  is 

a  vector  describing  the  translation  associated  with  this  symmetry  operation, 

At  this  point  in  the  program,  new  scattering  factor  coefficients 
are  read,  if  desired.  We  must  emphasize  two  points  of  caution  concerning 
the  changing  of  the  scattering  factors.   First,  since  the  blocking  of 
structure  factors  depends  upon  the  number  of  scattering  factors  used, 
if  that  number  is  changed,  new  structure  factors  must  be  read.  Violation 
of  this  rule  will  lead  to  unpredictable,  but  probably  catastrophic, 
results.  The  second  rule  requires  that  the  scattering  factors  be  changed 


only  when  file  ISC2  is  being  written,  that  is,  either  on  the  first  cycle 
of  a  given  run  or  when  new  structure  factors  are  being  read  from  cards. 
If  this  restriction  is  violated,  the  new  curves  are  loaded  to  core  and 
subsequently  included  in  the  data  set  for  this  structure,  but  they  are 
never  used  to  compute  scattering  factors  for  the  reflections . 

The  new  atom  parameters,  if  required,  are  now  read  from  cards. 
The  scale  factors  are  read  first,  followed  by  the  extinction  parameters  and 
the  atom  cards .  The  atom  parameters  needed  are  an  identification  label 
of  6  characters,  an  equipoint  fraction  which  is  the  multiplicity  of 
the  position  occupied  by  this  atom  divided  by  the  multiplicity  of  an 
atom  in  the  general  position  of  this  space  group,  the  total  occupancy 
of  the  site  (usually  one) ,  the  occupancy  of  species  A  (needed  if  more 
than  one  species  present  in  this  site) ,  the  identification  of  scattering 
factor (s)  to  be  applied  to  this  atom,  fractional  coordinates  for  the 
atom,  and  temperature  factor  information,  as  well  as  third  and  fourth 
cumulant  parameters,  consistent  with  the  temperature  factor  type  card 
described  below.  The  forms  of  the  cumulant  expressions  used  in  this 
program  are: 

F  =  Jf  expti^xJ  -  ZZIh^h^)  -  Zl(hkh^£  -WhhH^)] 
J  J  kJlm  k£  mn  J 

atoms 

•         f.i- 

x~  is  the  k fractional  coordinate  of  atom  j ,  3  is  the  array  of 

anisotropic  thermal  coefficients ,  and  c_  and  d_  are  the  arrays  of  third 
and  fourth  cumulant  coefficients .  After  the  atom  parameter  cards  have 
been  read,  the  special  parameter  cards  are  read.  The  first  card  gives 
the  total  number  of  special  parameters  (up  to  60) .  Then  one  card  for 
each  special  parameter  contains  its  name,  three  integers  specifying 
its  type,  and  its  initial  value. 

At  this  point  in  the  program,  the  parameter  refinement  selection 
information  is  read.  A  parameter  of  the  model  is  varied  in  the  least- 
squares  refinement  if  its  associated  indicator  is  one;  parameters  with 
zero  indicators  are  changed  only  if  they  are  linked  to  refined  parameters 
by  dependency  cards.  If  new  parameter  selection  information  is  read, 
constraint  data  are  also  read.  If  no  dependency  relationships  exist 
for  this  structure,  a  blank  card  must  be  included  in  the  data  deck. 
The  method  used  to  treat  parameter  constraints  is  from  Finger  (1969) 
and  essentially  consists  of  the  application  of  the  chain  rule  to  the 
derivatives  of  the  structure  factor  with  respect  to  the  parameters  of 
the  model.  The  user  must  determine  the  dependent  parameters  of  the 
model  using  the  rules  of  Levy  (1956)  for  thermal  coefficients  or  the 
equations  of  Finger  (1969)  for  chemical  constraints  and  then  prepare 
the  data  cards .  The  program  currently  will  not  allow  the  constraining 
of  the  total  occupancy  for  a  site  with  only  one  species  present;  however, 
this  may  be  accomplished  by  including  a  null  scattering  curve  which 
is  made  the  second  species  present.  The  amount  of  the  first  species 
may  then  be  refined  and  constrained. 


At  this  point,  all  structural  data  except  the  observed  structure 
factors  have  been  read.   If  these  are  to  be  read,  a  data  card  containing  a 
FORMAT  specification  is  read.  Two  optional  cards  may  be  read,  one  to 
specify  the  wavelength  if  it  is  needed  for  the  computation  of  extinction 
corrections  and  the  other  to  specify  an  "ignorance  factor",  p,  to  be  used 
in  computing-weights  according  to  the  formula 

w  =  l/[a2+  (pF   )2]. 
obs 

Subsequently,  the  subroutine  INPUT  is  called  to  read  and  process  a  reflec- 
tion card.   The  data  are  then  added  to  a  buffer  which,  if  full,  is  output 
on  file  ISC2 .   If  the  old  structure  factors  are  desired  and  this  is 
the  first  cycle  of  a  given  computer  run,  the  buffer  is  filled  from  file 
ISC1,  the  scattering  factors  are  recomputed  if  necessary  and  the  buffer 
is  written  on  file  ISC2.   If  this  is  not  the  first  cycle  of  a  given 
run,  the  old  structure  factors  are  read  from  file  ISC2 .  After  the  data 
for  a  reflection  have  been  obtained,  the  calculated  structure  factor 
and  partial  derivatives  are  computed  in  subroutine  SFAC,  the  results 
tested  to  determine  if  this  reflection  should  be  used  in  the  normal 
equations  matrix,  and  the  contributions  of  this  reflection  to  the  residuals 
are  calculated.  Next,  the  structure  factors  are  printed,  if  desired, 
with  the  rejected  data  marked  with  a  single  asterisk  in  the  right  margin. 
In  addition,  the  reflections  having  intensities  less  than  the  minimum 
observable  value  are  marked  with  a  single  asterisk  next  to  the  value 
of  F  ,   and  a  double  asterisk  in  the  right  margin. 

If  the  reflection  is  to  be  included  in  the  matrices,  subroutine 
MATRIX  is  entered;  then  the  next  observation  is  processed.  As  soon 
as  the  end  of  the  reflection  list  is  encountered,  the  residuals  for 
this  cycle  are  printed.  Next,  the  normal  equations  matrices  are  converted 
to  correlation  matrices  to  prevent  overflow  or  underflow  of  the  determinant 
during  the  matrix  inversion,  which  is  accomplished  by  subroutine  SYMINV, 
adapted  from  Busing  and  Levy  (1962).   If  the  matrix  is  non-singular, 
the  inverse  is  converted  to  a  variance-covariance  matrix  and  the  parameter 
shifts  are  computed.   Subroutine  RESET  is  entered  to  compute  the  new 
values  of  any  dependent  parameters.  The  ratios  of  the  parameter  shifts 
to  the  standard  deviations  are  tested  to  find  the  largest  absolute  value 
and  the  magnitudes  of  the  ratios  are  averaged.  At  this  point  in  the 
program,  the  correlation  matrix,  if  requested,  is  printed.   The  twenty 
off-diagonal  elements  of  the  correlation  matrix  with  largest  absolute 
values  are  always  printed.  Then  the  optional  bond  functions  are  computed 
by  routine  B0DAN1  and  the  thermal  functions  by  routine  ELVIB1. 

Finally,  a  run  continuation  card  is  read.   If  column  1  is  non-zero, 
another  cycle  is  to  be  performed.  However,  if  this  column  is  blank 
or  zero,  the  program  prepares  to  terminate.   If  the  new  values  for  the 
parameters  are  not  to  be  saved,  the  program  exits.  On  the  other  hand, 
if  these  parameters  are  to  be  kept,  a  new  label  record  is  written  on 
ISC1,  the  structure  factor  records  are  copied  from  ISC2  to  ISC1  and 
the  variance-covariance  record  is  output .  The  current  values  of  all 
parameters  may  be  punched  on  cards .  This  routine  then  returns  control 
to  the  main  program. 


Subroutine  B0DAN1 .   This  subroutine  uses  the  symmetry  operators 
supplied  on  cards  combined  with  the  action  of  centers  of  symmetry,  if 
present,  and  the  effects  of  multiple  lattice  points  and  cell  translations 
to  generate  all  atomic  locations  within  a  spherical  shell  centered  on 
each  atom  in  the  asymmetric  unit.  The  search  for  the  contents  of  this 
shell  is  done  efficiently  by  enclosing  the  central  atom  with  a  unit 
cell  shaped  parallelopiped  tangent  to  the  sphere.   For  each  location 
inside  this  box,  a  table  of  generated  positions  is  checked  to  determine 
if  this  location  is  unique  and  if  it  is,  the  bond  distance  is  computed. 
If  this  is  within  the  desired  range,  the  distance  and  coordinates  are 
printed  and  this  location  is  added  to  the  table  of  generated  positions. 

If  bond  angles  are  also  desired,  each  atom  in  the  asymmetric  unit 
is  used  as  the  central  atom  with  the  atom  positions  generated  in  the 
bond  distance  calculation  considered.   The  peripheral  atoms  are  identified 
only  by  their  alphanumeric  tag  and  a  number  which  points  to  an  entry 
in  the  bond  distance  output  where  the  coordinates  of  the  generated  atom 
are  listed.   Caution — if  the  maximum  bond  distance  is  large,  very  many 
angles  will  be  generated. 

Subroutine  ELVIB1 .   For  each  anisotropic  atom,  the  equivalent  iso- 
tropic temperature  factor  of  Hamilton  (1959)  ,  the  root-mean-square  (rms) 
amplitudes  of  vibration  (in  A)  of  the  principal  axes  and  the  orientations 
of  these  axes  with  respect  to  the  direct  cell  axes  are  computed,  or  if 
the  temperature  factor  matrix  of  an  atom  is  nonpositive-definite,  this 
is  reported.   This  routine  will  also  compute  the  thermal  corrections 
to  bonds  using  the  formulae  of  Busing  and  Levy  (1964) .   To  assist  the 
interpretation  of  the  orientations  of  thermal  ellipsoids  relative  to 
the  bonding,  the  rms  amplitude  of  the  ellipsoid  in  the  direction  of 
the  bond  is  computed  and  printed.   B0DAN1  is  used  to  generate  the  necessary 
atomic  locations . 

Subroutine  INPUT.   This  subroutine  reads  a  reflection  data  card, 
computes  sinz0/Az,  evaluates  the  scattering  factors  for  this  observation 
and  computes  a-p   f.or  the  observation.  Maximum  and  minimum  values  of 
Fobs  anc*  sin26/X2,  are  determined  in  this  routine. 

Subroutine  MATRIX.   This  subroutine  calculates  the  derivatives 
of  the  residual  from  the  derivatives  of  the  calculated  structure  factor 
and  accumulates  the  sums  required  in  the  normal  equations  matrices. 

Subroutine  MODIFY. 


and 


to  the  calculated  derivatives.  A  and  jB  are  the  real  and  imaginary  portions 
of  the  calculated  structure  factor,  P  is  the  independent  parameter,  p 


This 

routine  applies  the  equations 

9A  ' 
3*n 

=  Z 
m 

9A_!Zs 
9p   9p 
m   n 

9B  ' 
3P„ 

=  E 
m 

_3B  ^m 
9p    9p 
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represents  a  dependent  parameter,  the  summation  is  over  all  the  dependen- 
cies, the  unprimed  derivative  denotes  the  one  calculated  without  regard 
for  dependencies  and  the  primed  derivative  is  needed  to  form  the  least- 
squares  matrices. 

Subroutine  RCALC.   This  subroutine  accumulates  the  terms  needed  to 
calculate  the  discrepancy  factors  which  are  computed  for  all  reflections 
and  separately  for  those  data  used  in  the  refinement.  A  weighted  residual, 
R  which  is  computed  using  the  formula 


2 


Kw    Zw  (|Fo|  -  |Fd  )' 


2 


w   = 

Z  w  I  Fo  I 


for  refinement   on  J_  or 

R    =   Zw(|Fo|2   -    |FC|2)2 

w   ' ' ' ' — L- 

Zw  \?o\h 

for  refinement  on  J_2,  is  determined  for  each  class  of  data.   In  addition, 
the  conventional  discrepancy  index 

R=E{|  F°l  -  |FcH 
Z  |Fo| 

is  computed  for  all  reflections  in  each  class.   In  addition,  the  ranges  of 
variation  of  F  and  sin28/X2  are  each  divided  into  eight  subranges  and  the 
sums  needed  for  R  are  accumulated.   (Note:   the  maximum  and  minimum  values 
are  dynamically  determined  during  any  cycle  in  which  new  structure  factors 
are  read  and  the  distribution  of  residuals  for  such  a  cycle  will  be  mean- 
ingless.)  If  the  weighting  scheme  is  correct,  the  expression  (l/s)  Zw(F  -F  )2 
where  the  sum  is  taken  over  a  subset  of  the  full  data  set  containing  s 
reflections,  should  be  independent  of  the  composition  of  the  subset.   If 
we  divide  this  expression  by  the  constant  (1/t)  ZwF2 ,  where  t  is  the 
total  number  of  reflections  in  the  data  set  and  the  sum  is  taken  over  the 
entire  set,  and  take  the  square  root  of  the  quotient,  the  result  should 
also  be  independent  of  the  composition  of  the  subset.   This  number,  which 
can  be  considered  the  average  contribution  to  weighted  R  for  reflections 
in  the  subset  (Note  that  the  expression  is  equal  to  weighted  R  if  s=t.) 
is  printed  in  the  output  for  each  subset. 

Subroutine  RESET.   This  subroutine  is  entered  once  for  special  para- 
meters and  once  for  each  atom  to  compute  the  updated  values  of  all  depend- 
ent parameters.  Linear  relationships  and  rigid-body  parameters  are  com- 
puted directly.  Other  non-linear  systems  of  constraints  may  be  handled 
by  the  user  through  subroutine  SPVAL,  supplied  as  a  dummy. 

Subroutine  SFAC.   This  subroutine  computes  the  structure  factor  and 
partial  derivatives  with  the  symmetry  operators  applied  to  reciprocal  space; 
i.e.  the  Miller  indices  are  transformed  rather  than  the  atomic  coordinates. 
If  H1  represents  the  transpose  of  the  vector  of  Miller  indices  and  X'  . 


represents  the  jth  symmetry  transformed  coordinate  for  the  n_th  atom,  then 
the  argument  of  the  trigonometric  functions  needed  to  calculate  the  con- 
tribution of  this  atom  to  the  structure  factor  may  be  written  as  HTX'nj. 
Utilizing  the  expression  for  X'  from  page  3, 

HTX  =  HT(RX  +  T)  =  HTRX  +  HTT. 
The  final  form  is  the  method  used  in  this  subroutine. 

As  the  routine  is  entered  for  each  reflection,  the  matrix  products 
H  Rj  =  H!  and  HTTj  are  computed  and  saved.  Then  the  structure  factors 
and  partial  derivatives  are  accumulated.  If  3rd  and  4th  cumulant  terms 
are  absent,  the  structure  factor 


with 


and 


F=Eaf   ^  T  A?  A 
n  n  n  j  nj  nj 


P  .  =  exp  [2tt1  (H,TX  +  HTT.)] 
nj     r  n      j 


for  anisotropic  atoms  or 


T  .  =  exp  (-H'^  3  H!) 
nj     r     j   n  3 


T  .  =  exp  (-B  sin29A2) 
nj         n 


for  isotropic  atoms,  the  n  summation  is  over  the  atoms  in  the  asymmetric 
unit,  the  j_  summation  is  over  the  symmetry  operators,  a  is  an  occupancy 
term,  f  is  the  value  of  the  scattering  curve  for  the  n*h  atom,  i_  is  /-I 
and  $  *s  the  matrix  of  anisotropic  thermal  coefficients.   When  anomalous 
dispersion  is  included,  the  scattering  factor  becomes 


f  =  f   +  f '  +  if" 
n    o     n     n 
n 


with  f_Q  the  wave  length  independent  part ,  f^  the  real  part  of  the  anomalous 
dispersion  and  f^  the  imaginary  part.   Including  this  expression  in  the 
structure  factor  equation,  applying  Euler's  rule  and  writing  F_  in  the  form 
F  =  A  +  iB,  we  find  that 


(f   +  f ')  ?  T  .cos2tt(H!TX  +  HTT  ) 
on    n'  j  nj       j  n      j 


n   n 

-  a  f "  ?  T  4sin2TT(H!TX     +  HTT.) 
n  n  j     nj  j     n  j  J 


and 


=  Z  Ta   (f 
nL  n   or 

+  a  f"  ?  cos  2tt  (H/TX  +  HTT.) 
n  n  3         n  n      j  J 


+  f ' )  E  Sin  2tt(h:TX  +  HTT.) 
n   -         j  n      j 


The  required  partial  derivatives  are  found  by  differentiating  each  of  the 
above  equations  with  respect  to  the  components  of  Xn  and  (3.= .   After  all 
sums  have  been  accumulated,  the  subroutine  MODIFY  is  called  to  perform  any 
chainrule  operations  required.   If  this  is  to  function  properly,  all 
derivatives  must  be  calculated  although  the  refinement  indicators  would 
suggest  that  certain  derivatives  are  not  needed. 

Subroutine  SFAC  also  includes  an  entry  point  named  EXCALC  which  cal- 
culates the  extinction  factor  of  Zachariasen  (1967)  as  extended  by  Coppens 
and  Hamilton  (1970),  and  revised  by  Nelmes  and  Thornley (1973) .   Using 
this  scheme,  the  structure  factor  is  corrected  for  extinction  by  the 
relationship 


F*=F  y1/2 
c    c 

y  =  (1  +  B(9)Fc2  r*) 


-1/2 


-2  e4  2 (  +  cos426)   dA 


3(6)  = 


n2c4V2A  sin  26  (1  +  cos229)  dy 


or  x-rays 


X2         dA     - 

-t—  for  neutrons 


V2  A  sin  20    dy 


X_  is  the  wavelength,  V  is  the  cell  volume,  A  is  the  transmission  factor, 
u_  is  the  linear  absorption  coefficient,  and  e_,  m  and  c_  have  their  usual 
meanings.  The  program  allows  three  different  forms  for  the  extinction: 
a)  isotropic  with  r*  refined,  b)  Type  I  anisotropic  with  r*  =_/D'ZD)-1'2 
with  zj[j  refined  and  c)  Type  II  anisotropic  with  r*  =  (N'W  N)     with_ 
W.  .  refined.   If  p  is  a  unit  vector  parallel  to  the  incident  beam  and  s 
is  a  unit  vector  parallel  to  the  diffracted  beam,  D  is  a  unit  vector  par- 
allel' to  p  x  s  and  N  =  p  x  D.-  The  reference  coordinate  system  for  Z  and 
W  is  attached  to  the  crystal  and  has  its  x-axis  parallel  to  the  inci- 
dent beam  and  its  z_  axis  vertical  when  all  diff Tactometer  angles  are  0. 

Subroutine  SYMINV.   The  method  used  by  this  matrix  inversion  routine 
has  been  described  by  Busing  and  Levy  (1962)  and  was  coded  by  Mr.  James 
Mundstock  of  the  Numerical  Analysis  Center  of  the  University  of  Minnesota. 

Subroutine  REJECT .   This  routine  may  be  used  to  exclude  reflections 
from  the  normal  equations  matrices  for  other  than  the  standard  reasons. 
If  a  reflection  is  to  be  rejected,  the  variable  NREJ  in  Common  Block  /F/ 
should  be  set  to  one. 
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Subroutine  WEIGHT .   This  routine  may  be  used  to  compute  a  non-standard 
weighting  scheme.   Upon  entry,  the  information  contained  in  the  standard 
deviation  field  of  the  structure  factor  card  will  be  in  the  variable  WT 
in  Common  Block  /F/.   This  variable  must  be  set  to  the  desired  weight 
information  when  the  routine  is  exited. 

Subroutine  RDEGN.   This  subroutine  reads  the  dependent  parameter  in- 
formation from  cards,  and  rearranges  the  array  of  relationships  in  ascend- 
ing order  first  of  independent  atom  and  parameter  designators  and  then  of 
dependent  atom  and  parameter  designators.   Special  parameters  are  arbitra- 
ily  assigned  atom  number  0. 

Subroutine  DERIVS .   This,  subroutine  is  entered  once  at  the  beginning 
of  each  cycle .  It  computes  the  current  values  of  derivatives  of  conven- 
tional parameters  with  respect  to  rigid  body  parameters.   (Other  non-linear 
derivatives  may  be  supplied  by  the  user  through  subroutine  SPDERI) .   If 
any  independent  parameters  are  themselves  actually  dependent  on  other  par- 
ameters (such  as  when  the  values  of  the  elements  of  the  rigid-body  motion 
tensors  are  restricted  by  symmetry)  the  derivatives  with  respect. to  the 
true  independent  parameters  are  computed.  Then  the  independent  parameter 
designators  are  replaced  by  their  row  indices  in  the  normal  equations 
matrix,  and  the  array  is  packed  to  eliminate  any  redundant  relations. 

Subroutine  PAREXT.  This  subroutine  reads  the  special  parameter 
cards,  and  prints  the  names  and  initial  values  of  the  parameters. 

Subroutine  TRNFRM.   This  subroutine  is  entered  once  at  the  beginning 
of  a  run  to  compute  the  transformation  matrices  required  to  convert  para- 
meters from  the  crystal  system  to  a  standard  orthonormal  coordinate  system 
and  back  again. 

Subroutine  FILL IN.  This  subroutine  stores  special  parameter  values 
in  the  proper  positions  in  rigid-body  parameter  arrays . 

Function  ERRCAL.   This  function  subroutine  determines  the  standard 
deviation  of  a  dependent  parameter  according  to  the  approximate  formula 

*    ,  n   n  9P,  3 P.     .1/2 
K    1=1  j=l  3P4  3P.   1J 

where  n  is  the  number  of  refined  parameters  and  cr-j-j  is  the  covariance  of 
the  ith  refined  parameter  with  the  j  th  refined  parameter. 

Subroutine  EULER.   This  subroutine  computes  an  orthogonal  matrix 
which  is  necessary  to  transform  parameters  from  a  general  orthonormal 
coordinate  system  to  the  standard  orthonormal  coordinate  system  defined 
by  x  parallel  to  a  and  z  parallel  to  c.*.   The  relative  orientations  are 
defined  by  Eulerian  angles  <j>,x>  and  u.   Here  oi  is  the  angle  through 
which  the  special  coordinate  system  must  be  rotated  clockwise,  as  viewed 
down  the  positive  z  axis  of  the  standard  system,  about  the  z  axis  of 
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the  standard  system  to  bring  the  z  axis  of  the  special  system  into  the 
x-z  plane  of  the  standard  system;  x  is  a  clockwise  rotation  about  the  y 
axis  of  the  standard  system  (viewed  down  the  positive  y  axis)  to  bring 
the  z  axes  of  the  two  systems  into  coincidence;  <J>  is  a  clockwise  rotation 
about  the  common  z  axis  to  bring  the  x  and  y  axes  into  coincidence. 
Three  additional  entry  points,  DTDPHI,  DTDCHI,  and  DTDOMG  compute  the 
derivatives  of  the  elements  of  the  transformation  with  respect  to  (j>,  x> 
and  a)  respectively. 

Subroutine  MULMM.   This  subroutine  forms  the  product  of  two  3x3 
matrices.  An  alternate  entry,  MULMV,  multiplies  a  three  component 
column  vector  by  a  3  x  3  matrix. 

Subroutine  TFM.   This  subroutine  has  two  entry  points.   TFM2  trans- 
forms a  second  cumulant  tensor  in  an  orthonormal  system  to  a  set  of  aniso- 
tropic temperature  factors  in  the  crystal  system.   TFM3  does  the  same  for 
a  third  cumulant  tensor. 

Subroutine  TRNSPS.   This  subroutine  generates  the  transpose  of  a 
3x3  matrix. 

Function  IPAK.    This  function  packs  4  integer  variables  into  a 
single  computer  word.  The  version  given  is  appropriate  only  to  the  UNIVAC 
1108.   Other  versions  must  be  used  for  other  machines.   If  the  integer 
word  length  of  the  machine  is  shorter  than  32  bits  other  adjustments 
must  be  made. 

Functions  SKTLS,  D2KDT,  D2KDL,  and  D2KDS.   These  functions  compute 
the  second  cumulants  for  an  atom  and  their  derivatives  with  respect  to  the 
rigid  body  motion  parameters  T,  L,  and  S/. 

Functions  TKTLS ,  D3KDL,  and  D3KDS.   compute  the  third  cumulants 
and  their  derivatives  with  respect  to  L  and  S. 

Subroutine  CORECT.   This  subroutine  applies  libration  corrections  to 
the  atomic  positions  in  a  rigid  body. 

Function  DETERM  and  subroutine  SOLVE  are  used  by  CORECT.   Subroutine 
VXV  forms  the  vector  product  of  two  vectors .   Functions  AIJ,  BIJK,  CIJKL, 
DIJKLM,  EIJ,  FIJ,  AXA,  AXAXB,  PIJKLM,  EL4 .  are  a  set  of  functions  used 
in  the  rigid  body  relations. 

Common  Blocks  and  Variables 

1.   Common  Block  /LSMAT/ 

B(200)  Normal  equations  vector 

IEXT  Extinction  type  indicator 

IPARA  (1000)  Atom  Parameter  refinement  indicators 

MRANK  Order  of  the  normal  equations  matrices 

KCYL  Cycle  number  for  the  data  set 

LABEL  (3)  12  character  label  for  the  data  set 
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JSCRF  (10) 

NRSF 

SIGMA 

MDEGN 

A(22740) 


Common  Block  /A/ : 

BETA(6,60) 

G(3,3) 

IS0T(60) 

NATOM 

TAG (2, 60) 

TITLE (18) 

XYZ(3,60) 

PI 

CR(10,60) 

DR(15,60) 


Refinement  switches  for  scale  factors 
Number  of  scale  factors  refined 
Standard  deviation  of  a  parameter 
Twice  the  number  of  constraint  relationships 
Constraints  array  and  normal  equations  matrix 
or  inverse.   The  dimensions  are  a  maximum. 
The  actual  number  for  a  given  problem  is 
given  by  the  expression  MDEGN+ (MRANK* (MRANK 
+l))/2.   If,  as  on  the  UNIVAC  1108,  the 
operating  system  allows  dynamic  core  alloca- 
tion, this  number  can  be  determined  at  run 
time . 


Temperature  factor  coefficients 
Real  space  metric  tensor 
Temperature  factor  type  indicators 
Number  of  atoms  in  the  asymmetric  unit 
Alphanumeric  atom  tags  or  labels 
Alphanumeric  title  information 
Atomic  positions 
Value  of  tt 

Third  cumulant  coefficients 
Fourth  cumulant  coefficients 


3.  Common -Block  /C/: 

CELLP(3,4) 

ICENT 

JCELL 

NSYM 

RMAT(3,3,24) 

TRANS (3, 24) 

4.  Common  Block  /D/i 


Coordinates  of  lattice  points 
Centric  indicator 
Number  of  lattice  points 
Number  of  symmetry  operators 
Rotational  symmetry  operators 
Translational  symmetry  operators 


ACALC 

AIND(60) 

ANOM(2,20) 

BCALC 

BIND (60) 

DADB(6,60) 

DADO (60) 

DADX(3,60) 
DBDB(6,60) 

DBDO(60) 
DBDX(3,60) 


Real  part  of  calculated  structure  factor 
Individual  atom  contributions  to  ACALC 
Anomalous  dispersion  corrections 
Imaginary  part  of  calculated  structure  factor 
Individual  atom  contributions  to  BCALC 
Derivatives  of  ACALC  with  respect  to  temper- 
ature factors 
Derivatives  of  ACALC  with  respect  to  occu- 
pancy 
Derivatives  of  ACALC  with  respect  to  position 
Derivatives  of  BCALC  with  respect  to  temper- 
ature factor 
Derivatives  of  BCALC  with  respect  to  occupancy 
Derivatives  of  BCALC  with  respect  to  position 
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H(3) 
IH(3) 

ISCAT(2,60) 

ISCL 

MFSEL 

OCCA(60) 

OCCUP(60) 

SCALE (10) 

SCAT (20) 

SITE(60) 
DADC(10,60) 

DBDC(10,60) 

DADD(15,60) 

DBDD(15,60) 

AMINUS 
BMINUS 
DRDE(6) 

NDEGN 

5.  Common  Block  /E/ : 

DMAX 
DMIN 

RECELL(6) 
EXANG(3) 

6.  Common  Block  /F/ i 

DELTA 

FCALC 

FOBS 

FOMAX 

FOMIN 

IREJ 

MREF 

MREJ 

NREJ 

NUMBER(36) 

RDEN(36) 

RHO 

RHOMN 

RHOMX 

RNUM(36) 


Miller  indices  in  real  representation 

Miller  indices  in  integer  representation 

Table  used  to  select  atomic  scattering  factor 

Number  of  scale  factor  for  this  reflection 

Function  selector 

Occupancy  of  species  A 

Total  occupancy  of  site 

Scale  factors 

Values  of  scattering  factor  for  current 

reflection 
Equipoint  'fraction  for  atom  site 
Derivatives  of  ACALC  with  respect  to  3rd 

cumulant 
Derivatives  of  BCALC  with  respect  to  3rd 

cumulant 
Derivatives  of  ACALC  with  respect  to  4th 

cumulant 
Derivatives  of  BCALC  with  respect  to  4th 

cumulant 
ACALC  for  opposite  absolute  orientation 
BCALC  for  opposite  absolute  orientation 
Derivatives  of  r*  (extinction  factor)  with 

respect  to  coefficients  of  tensor 
Number  of  constraint  relationships 


Maximum  bond  distance 
Minimum  bond  distance 
Reciprocal  cell  lengths  and  cosines 
4>,  X  and0+  a)  (in  radians)  for  current 
reflection 


Weighted  residual  for  current  reflection 

Calculated  structure  factor 

Observed  structure  factor 

Maximum  observed  structure  factor 

Minimum  observed  structure  factor 

Rejection  mode  indicator 

Refinement  mode  indicator 

Minimum  observable  indicator 

Used  to  indicate  whether  or  not  a  reflection 

is  rejected 
Counters  for  subdivisions  of  R  factors 
Denominators  for  R  factors 
Current  value  of  sin20/X2 
Minimum  value  of  RHO 
Maximum  value  of  RHO 
Numerators  for  R  factors 
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WT 

EXBETA 

EXX 

EXY 

FISGN 

IFSGN 

EXTIN(6) 

KEXTRF(60) 

DELMAX 

7.  Common  Block  /G/ : 

BUFF (1023) 

FCURVE(9,20) 
FRMT(18) 

IWT 

NSCAT 

BETFAC 

WVLNTH 

WTFAC 

8.  Common  Block  /TAPE/: 


Weight  for  current  reflection 

Factor  (3(9)  for  extinction  correction 

Factor  B(9)F2r*  used  in  extinction 

-1/5 
(1  +  EXX)   ' 

Absolute  orientation  parameter 

Orientation  determination  indicator 

Extinction  parameter  matrix 

Refinement  switches  for  extinction  parameters 

Maximum  permitted  value  of  AF/ct 


Array  used  primarily  for  blocking  structure 

factor  records 
Scattering  curve  coefficients 
Variable  format  specification  for  reading 

structure  factors 
Weight  calculation  ,type  indicator 
Number  of  scattering  factors 
Used  in  extinction  correction 
Wavelength  of  radiation 
"Ignorance  factor"  in  computation  of  weights 


IN 

INN 

IOUT 
ISC1 

ISC2 

ISC3 
ISC4 
ISC5 

IPUN 


Logical  unit 
Logical  unit 

input  unit 
Logical  unit 
Logical  unit 

file  (10) 
Logical  unit 

file  (29) 
Logical  unit 
Logical  unit 
Logical  unit 

list  (28) 
Logical  unit 


number  of  the  card  read  unit  (5) 
number  of  the  structure  factor 
(5  or  11) 

number  of  the  print  file  (6) 
number  of  the  least-squares 

number  of  the  intermediate 

number  of  the  auxiliary  file  (11) 
number  of  the  Fourier  file  (12) 
number  for  saving  constraints 

number  for  punch  file  (1) 


Dependent  Parameter  Input 

In  certain  structures,  the  parameters  of  the  model  may  not  be  inde- 
pendently variable.   If  they  are  merely  constrained  to  have  a  fixed  value, 
there  are  no  difficulties  as  the  refinement  selection  variable  for  this 
parameter  is  set  to  zero.   However,  if  a  parameter  of  the  model  is 
required  to  be  equal  to  a  linear  combination  of  variable  quantities,  the 
situation  is  greatly  complicated.   In  the  least-squares  refinement  of 
such  a  structure,  the  contribution  of  the  dependent  parameter  to  the 
derivative  of  the  structure  factor  with  respect  to  the  independent 
variable  must  be  computed  if  the  shifts  are  to  be  properly  evaluated. 
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Similarly,  the  correct  value  of  the  dependent  parameter  must  be  restored 
at  the  end  of  the  cycle  if  the  structure  factors  are  to  be  properly 
computed  in  the  next  cycle .   This  program  utilizes  a  new  algorithm  to 
perform  the  derivative  modification  and  parameter  restoration.  The  method 
of  preparing  the  input  data  cards  will  be  shown  with  the  aid  of  examples . 

Example  1:  Atom  6  is  an  anisotropic  atom  on  a  3-fold  axis  parallel 
to  c_.  "  Employing  the  rules  of  Levy  (1956) ,  the  following  relationships 
are  required  for  this  atom: 


13   M23 


=  0 


22   2  312 


The  first  equations  relate  parameters  to  a  constant  and  thus  may  be 
ignored.   However,  this  is  not  the  case  for  the  second,  and  since  3n 
must  be  the  independent  variable,  the  following  equations  are  applied  as 
constraints  to  the  least-squares  solutions: 


22 
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312  =  1/2  3n. 

The  resulting  Parameter  Constraint  Cards  (See  page  24)  for  this  situation 
would  be: 

Dependent    Independent    Constant    Comments  (not  on  card) 

6   6         6   5        1.0 


0.5 


622  =  3-1 1  for  atom  6 

312  =  o.5  eu 

End  of  list 


Example  2:   Grunerite  has  four  independent  sites  (atoms  1  to  4)  for 
Fe  and  Mg,  but  the  total  amount  of  each  is  known  from  the  chemical  analy- 
sis.  Therefore,  we  want  to  refine  the  values  of  the  Mg-occupancies  for 
the  first  three  atoms  and  constrain  the  fourth  to  agree  with  the  bulk 
chemistry.   In  the  notation  of  Finger  (1969) ,  the  constraining  equation  is 

E  b  a   =  C 
m  m  mn    n 

with  bju  the  multiplicity  of  the  mth  site,  ajnn  the  fractional  occupancy 
for  the"  nth  species  in  the  mth  site,  and  C  ""the  total  number  of  atoms  of 
species  n  per  unit  cell.   In  grunerite,  sites  1  through  4  have  multiplici- 
ties 4,  4,  2,  and  4  respectively  and  there  are  two  formula  units  per  unit 
cell.   If  we  constrain  the  formula  to  be  Feg  2^0.8^8^22 ^*^2  anc^  su^~ 
stitute  into  the  above  equation,  we  find  that 
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Solving  for  a 


4Mg» 


4alMg  +  4A2Mg  +  2a3Mg  +  4a4Mg  =  1*6 


a4Mg  =  °'4  "  alMg"  a2Mg  "  1/2  a3Mg' 


The  Parameter  Dependency  Cards  needed  to  constrain  this  chemistry  would  be: 

Comments  (not  on  card) 

Initialization 

Site  1  dependence 

Site  2 

Site  3 

End  of  list 


Dependent 

Indep< 

sndent 

Constant 

4        1 

0 

0 

0.40 

4        1 

1 

1 

-1.0 

4        1 

2 

1 

-1.0 

4        1 

3 

1 

-0.5 

0 

Alternatively,  the  latest  version  allows  up  to  5  independent  parameters 
to  affect  one  dependent  one  on  the  same  card,  so  the  first  four  cards 
above  may  be  replaced  by  the  card  as  follows : 

4  10  0  0.40  1  1  -1.0  2  1  -1.0  3  1  -0.5 

Note  that  when  occupancies  are  constrained,  the  dependent  value  must 
first  be  initialized  with  a  card  containing  zeros  in  the  independent  par- 
ameter fields. 

Example  3 :  Like  many  chemical  groups ,  phosphate  groups  often  have 
thermal  vibration  that  can  be  represented  by  the  motions  of  a  rigid  body, 
which  can  be  expressed  as  the  elements  of  the  T,  L,,  and  S_  tensors  intro- 
duced by  Schomaker  and  Trueblood  (1969) .   Expressions  for  the  second  and 
third  cumulants  of  the  probability  distribution  functions  for  individual 
atoms  (Johnson,  1969)  as  functions  of  the  elements  of  £  and  L   are  given 
by  Prince  and  Finger  (1973)  ,  and  the  corresponding  expressions  including 
_S  can  be  readily  derived  (Prince,  unpublished  but  given  as  a  FORTRAN 
expression  in  this  program) .   If  atom  3  is  an  oxygen  atom  of  a  phosphate 
group  and  special  parameters  1  through  24  are  3  coordinates  of  the 
origin  to  which  the  motions  are  referred,  6  elements  of  ^T,  6  elements  of 
L,  and  9  elements  of  S,  the  constraint  cards  will  be  set  up  as  follows: 
First,  because  of  the  so  called  "trace  of  S  singularity",  the  three 
diagonal  elements  of  £  are  not  independent,  and  S33  can  be  set  equal  to 
"Sll  "  S22  by  means  o?  the  card  0  18    0  16    -1.    0  17    -1. 
The  anisotropic  temperature  factors  may  depend  on  all  of  the  rigid 
body  elements,  so  there  must  be  six  cards  of  the  form 


3     5 

0 

4 

0. 

0-24 

0 

3     6 

0 

4 

0. 

0-24 

0 

3     7 

0 

4 

0. 

0-24 

0 

3     8 

0 

4 

0. 

0-24 

0 

3     9 

0 

4 

0. 

0-24 

0 

3  10 

0 

4 

0. 

0-24 

0 
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Here  the  minus  sign  before  the  24  indicates  that  the  dependent  parameter 
depends  on  special  parameters  4  through  24  inclusive.   The  coefficients 
are  given  as  zeros  because,  being  non-linear,  they  change  from  cycle  to 
cycle  and  are  therefore  computed  in  each  cycle  according  to  the  current 
values  of  the  refined  parameters.   Third  cumulants  are  independent  of  T, 
but  may  be  dependent  on  any  element  of  L  or  S.   Therefore  a  typical 
constraint  card  for  a  third  cumulant  would  have  the  form 

3  12    0  10    -0.    0-24    0. 

The  program  uses  the  information  on  the  constraint  cards  to  allocate 
space  in  any  array  and  then,  in  non-linear  cases,  fills  in  the  space  at  the 
beginning  of  each  cycle.   Therefore  if  a  given  relationship  is  not  listed 
it  will  never  be  used  whether  or  not  its  value  is  actually  zero.   Because 
of  symmetry  some  relationships,  particularly  those  involving  T,  may  be 
identically  zero.   It  will  save  space  and  almost  certainly  time  to  omit 
these  from  the  list. 

Details  of  Data  Input 

1.  File  Name  Card:   Columns  1-12  contain  an  alphanumeric  label  to  be 

written  into  a  new  data  set.  This  label  will  be 
printed  when  the  file  is  read . 

2.  Title  Card:  Any  desired  Hollerith  information  in  columns  1-72.   This 

information  appears  at  the  top  of  most  pages  of  the  out- 
put. 

3.  Control  Card: 

Column   Variable  Description 

1  MFSEL       ' Function  selector 

0  -  Structure  factors  only 

1  -  Least-squares  refinement 

2  MSF         Structure  factor  output  selector 

0  -  No  structure,  factor  output 

1  -  Printed  structure  factors 

2  -  Printed  and  written  (for  a  Fourier  program) 

on  file  ISC4 

3  -  Written  on  file  ISC4 

3  MCORR       Correlation  matrix* 

0  -  Not  printed 

1  -  Printed 

4  MBODAN       Bond  distance,  angle  selector 

0  -  None  computed 

1  -  Distances  only 

2  -  Distances  and  angles  computed 

5  MELVIB       Ellipsoids  of  vibration 

0  -  None  computed 


*The  values  of  the  20  largest  off-diagonal  correlation  coefficients  will 
be  printed  for  every  least-squares  cycle. 
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3.   Control  Card  (Cont.) 


Column   Variable 


Description 


8 

NEW(l) 

9 

NEW(2) 

10 

NEW(3) 

11 

NEW  (4) 

12 

NEW  (5) 

13 
14 

15 


1  -  Ellipsoids  computed  for  anisotropic  atoms 

2  -  Ellipsoids  and'  thermal  corrections  to  bonds 

computed 

6  MREF        Refinement  type  indicator 

0  -  Refine  on  F„ 

1  -  Refine  on  F 

7  IFILE        Rewrite  data  set  indicator 

0  -  Data  set  rewritten  if  this  is  the  last  cycle 

of  the  current  run 

1  -  Data  set  not  rewritten 

2  -  Data  set  rewritten.  On  input  the  standard 

deviations  of  parameters  refined  in  previous 
cycles  will  be  read  in. 

3  -  Data  set  not  rewritten.  ,  Standard  deviations 

read  in. 
If  non-zero ,  read  new  cell  parameters* 
If  non-zero,  read  new  bond  distance  limits 
If  non-zero,  read  new  symmetry  operators 
If  non-zero,  read  new  scattering  factorst 
If  non-zero,  read  new  scale  factors  and  atom  para- 
meters.  If  equal  to  one,  all  atoms  will  be  read 
If  equal  to  two,  only  selected  atoms  will  be  read. 
If  equal  to  3  only  scale  and  extinction  will  be 
read.   If  equal  to  4  only  special  parameters  will 
be  read. 
NEW(6)       If  non-zero,  read  temperature  factor  type  card 
NEW(7)       If  non-zero,  read  parameter  selection  card(s)  and 

parameter  dependency  card(s) 
NEW(8)       If  non-zero,  read  new  structure  factors.  Values  1, 
2 ,  and  3  correspond  to  various  options  for  the 
order  of  the  list,  as  explained  below. 


For  all  zero  values  in  columns  8-15 ,  the  old  values  are  read  from  file 
ISC1  or  the  values  in  memory  are  used .  The  above  integers  are  required  for 
all  cycles. 


16-20    DMIN 
21-25    DMAX 


Minimum  bond  distance  in  Angstroms.  Needed  only  if 

NEW(2)  4   0. 
Maximum  bond  distance  in  Angstroms.  Needed  only 

if  NEW(2)  4   0. 


*Note  that  if  the  new  cell  parameters  will  affect  sinS/X,  and  therefore 
the  atomic  scattering  factors  for  various  reflections,  structure  factors 
must  be  read  in. 

tWarning:  New  scattering  factor  curves  must  be  read  only  in  the  first 
cycle  of  a  given  computer  run  or  when  new  structure  factors  are  being 
read  in,  or  they  will  never  be  used  in  computation. 
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3.   Control  Card  (Cont.): 

Column   Variable  Description 

26  ICENT       Centric  indicator,  needed  only  if  NEW(3)  4   0 

0  -  Centric 

1  -  Acentric 

27  ICELL        Cell  type  indicator,  needed  only  if  NEW(3)  4   0 

0  -  Primitive  cell 

1  -  A-centered  cell 

2  -  B-centered  cell 

3  -  C-centered  cell 

4  -  I-centered  cell 

5  -  F-centered  cell 

6  -  R-centered  cell 

28  MANOD       Anomalous  dispersion  indicator,  needed  only  if 

NEW(4)  4   0 

0  -  No  corrections  are  to  be  used 

1  -  Include  anomalous  dispersion  corrections 
29-30     NSCAT       Number  of  scattering  curves.   This  is  needed  only 

if  NEW (4)  ^0.   (0  <  NSCAT  <  21) 
31-32     NATOM       Number  of  atoms  in  asymmetric  unit .   This  is  needed 

only  if  NEW (5)  ^  0.   (0  <  NATOM  <  61) 
33-34     NSF        Number  of  scale  factors.   This  is  needed  only  if 

NEW (5)  4   0.   (0  <  NSF  <  11) 

35  NEWEXT       Extinction  change  indicator  (needed  only  if  NEW(5) 

4   0) 

0  -  New  extinction  parameters  not  read 

1  -  Extinction  parameters  read. from  cards 

36  IWT         Weight  type,  needed  only  if  NEW(8)  4   0. 

0  -  Standard  deviation  of  F  on  cards 

1  -  Standard  deviation  of  F2  on  cards 

2  -  Weight  on  cards 

3  -  Sigma  for  F  is  to  be  computed  in  WEIGHT 

4  -  Compute  weight  in  subroutine  WEIGHT 

5  -  Standard  deviation  of  F  on  cards  will  be 

altered  by  "ignorance  factor"  by 
a     =  (a2  +(PF)2)!/2 

37  INN         Structure  factor  source  indicator,  needed  if 

NEW(8)  4   0. 

0  -  Structure  factors  on  cards 

1  -  New  structure  factors  from  file  ISC3 

38  IREJ        Rejection  parameter,  needed  only  if  NEW(8)  4   0. 

0,5  -  No  reflections  rejected 

1.6  -  Reject  if  I  less  than  the  minimum  obs. 

2.7  -  Reject  if  I  less  than  the  minimum  obs. 

or  |  Fo-Fc  |  hy   greater  than  DELMAX 

3.8  -  Reject  if  |Fo-Fc|/a  greater  than  DELMAX 

In  both  of  these  cases  weights  will  be 
a^red  by  w'  =  wCl^x^x^)  ,  where 
x  =  | Fo-Fc | / (aF*DELMAX) 

4.9  -  Use  special  routine  REJECT 
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3.   Control  Card  (Cont.) 


Column   Variable 


39-43 


61 


NOTE: 


DELMAX 


44-53 

JSCRF 

54-59 

KEXTRF 

60 

JXPAR 

IFSGN 


Description 

Values  of  IREJ  from  5  to  9  allow  the  rejection  mode 
to  be  changed  even  though  NEW(8)  =  0. 

Maximum  allowable  difference  between  observed  and 
calculated  structure  factors  on  an  absolute 
scale.  This  is  needed  only  if  NEW (8)  4   0  and 
IREJ  is  2  or  3,  or  if  IREJ  is  7  or  8. 

Scale  factor  refinement  switches  (needed  only  if 
NEW(7)  4   0) 

Extinction  factor  refinement  switches  (needed  only 
if  NEW(7)  4   0) 

Extra  parameter  indicator 

0  -  No  extra  parameters 

1  -  Extra  parameters  will  be  read 

(needed  only  if  NEW(5)  4   0) 
Orientation  determination  parameter 

0  -  Do  not  determine  the  orientation 

1  -  Change  to  orientation  with  lower  weighted 

r  (needed  only  if  NEW(7)  4   0) 


4.   Cell 

Card :  Ne< 

Columns 

Variable 

1-8 

A(l) 

9-16 

A(2) 

17-24 

A(3) 

25-32 

A(4) 

33-40 

A(5) 

41-48 

A(6) 

Needed  only  if  NEW(l)  4   0. 


Description 


a_  axis  length  in  angstroms 

b_  axis  length 

c_  axis  length 

Alpha  in  degrees 

Beta  in  degrees 

Gamma  in  degrees 


5.   Symmetry  Cards:   Needed  only  if  NEW(3)  4   0, 


Columns 

Variable 

1-3 

JTNS (1) 

5-6 

IRMAT(l.l) 

8-9 

IRMAT(2,1) 

11-13 

JTNS (2) 

15-16 

IRMAT(1,2) 

18-19 

IRMAT(2,2) 

21-23 

JTNS (3) 

25-26 

IRMAT(1,3) 

28-29 

IRMAT(2,3) 

34 

IEF 

Description 

Translational  part  of  x_' 

First  positional  part  of  x' 

Second  positional  part  of  x.' 

Translational  part  of  v_' 

First  positional  part  of  y_' 

Second  positional  part  of  y_' 

Translational  part  of  z/ 

First  positional  part  of  z/ 

Second  positional  part  of  z_' 

End  flag,  non-zero  on  last  card  only 

0  -  Another  symmetry  card  follows 

1  -  This  is  the  last  card 


20 


The  translational  operators  are  of  the  following  form: 

bbb  1/2  1/3  2/3  1/4  3/4  1/6  or  5/6.  The  positional  operators  are  of 
the  form  bW  =  W  or  -W  with  W  equal  to  X,  Y,  or  Z  and  b  signifies  a 
blank.  The  symmetry  position  X  Y  Z  must  be  included. 

6.   Scattering  Factor  Cards:   Needed  only  if  NEW(4)  4   0. 


These  coefficients  must 
satisfy  the  equation 

*  2 

f (x)  =  c  +  E  a  exp(-b  x  ) 

i=l 


Columns 

Description 

1-8 

c 

9-16 

17-24 

25-32 

33-40 

41-48 

49-56 

57-64 

65-72 

73-80 

Identifier 

with  x  =  sine A 


Coefficients  for  scattering  curves  of  this  form  have  been  tabulated 
by  Cromer  and  Waber  (1965)  for  relativistic  Dirac-Slater  wave  functions, 
by  Cromer  and  Mann  (1968)  for  Hartree-Fock  wave  functions  and  by  Doyle 
and  Turner  (1968)  for  relativistic  Hartree-Fock  wave  functions.   If 
anomalous  dispersion  coefficients  are  being  used  (MANOD  =  1) ,  the  above 
card  is  followed  by  one  with  the  real  part  of  the  anomalous  dispersion 
in  col.  1-8  and  the  imaginary  part  in  col.  9-16.  Coefficients  for  the 
hydrogen  atom  are  given  by  Forsyth  and  Wells  (1959) . 

NOTE:  The  scattering  curves  should  be  changed  only  on  the  first 
cycle  of  a  given  computer  run  or  when  NEW(8)  =  1.  New  structure  factors 
must  be  read  (NEW(8)  =  l)if  the  number  of  scattering  curves  is  to  be 
changed . 

7.  Scale  Factor  Cards:  Needed  only  if  NEW(5)  4   0. 

The  NSF  scale  factors  are  punched  with  8  columns  per  value  and  9 
values  per  card  and  should  be  the  values  needed  to  place  the  calculated 
structure  factors  on  the  same  scale  as  the  observed . 

8.  Extinction  Cards:  Needed  only  if  NEW (5)  4   0  and  NEWEXT  4   0.   Format 
(6E12.6) 

Columns     Value 

1-12  r*  for  isotropic  extinction  or  r*  for  anisotropic 

13-24  r* 

25-36  r*it 

37-48  r*f;J 

49-60  r*r^ 

61-72  r*g 
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9.   Atom  Parameter  Cards:  Needed  only  if  NEW (5)  4   0. 

If  NEW (5)  =  1,  the  entire  atom  list  is  read  with'  two  cards  per  atom. 
If  NEW(5)  =  2,  the  Atom  Number  is  read  with  Format  (12)  followed  by 
the  two  cards  for  that  atom.   Card  reading  is  terminated  by  a  blank 
Atom  Number  card. 


Card 


Columns 

1-6 

7-12 

13-18 

19-24 


Card  2  1-2 
3-4 

9-15 
16-22 
23-29 
30-36 

37-43 

44-50 
51-57 
58-64 
65-71 
73 


Variable 

TAG 
SITE 

occur 

OCCA 

ISCAT(l) 
I SCAT (2) 

XYZ(l) 
XYZ(2) 
XYZ(3) 
BETA(l) 

BETA (2) 

BETA (3) 
BETA (4) 
BETA (5) 
BETA (6) 
IEF 


1-8 
9-16 
17-24 
25-32 
33-40 
41-48 
49-56 
57-64 
65-72 
73-78 
80 


:iu 

:222 
:333 

;m 
;i22 
;ii3 
;i33 

:223 
:233 


Description 

Alphanumeric  identifier 
Equipoint  fraction 
Total  occupancy  of  site 
Occupancy  of  species  A  in  this  site 

Needed  only  if  multiple  occupancy 
Number  of  scattering  curve  for  species  A 
Number  of  scattering  curve  for  species  B 

(blank  if  single  element  in  site) . 
x  of  atom  in  fractional  coordinates 

Z 
z_ 

B_  for  an  isotropic  atom  or  (3, ,  for  an  aniso- 
tropic atom 

622 

e33 

312 

Jl3 

P23 

Indicator  of  third  cumulant 


Card  3  -  Needed  if  IEF  4   0, 


Fourth  cumulant  ind. 


Card  4-5  -  Needed  if  IEF  4   0, 


Punch  dim,  d2222,  d3333,  dm2,  dm3,  dm2,  dn33,  d1222>  d1223> 
d  2   ,  d..-.,  d222  ,  d 2233»  d2333  in  order  with  8  columns  per  value. 
10.   Special  parameter  cards:   Needed  only  if  NEW(5)  4   0  and  JXPAR  4   0. 

Card  1:   Columns  3-4  contain  the  number  of  special  parameters,  NPAR. 
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If  NPAR>0  the  succeeding  cards  contain  the  names,  initial  values,  and 
identifiers  according  to  the  following  form. 


Columns 
Columns 


1-8 
9-10. 


Columns  11-12 


Column  13-14 


Name  of  parameter,  up  to  8  alphanumeric  characters. 
Type  of  parameter,  0- linear  function  with  coefficient 

given  on  constraint  card. 
1-rigid  body  system 

2-defined  by  user  in  subroutines  SPVAL  and  SPDERI 
Integer  indicating  which  rigid-body  system,  if  more  than 

one ,  this  parameter  refers  to . 
Parameter  designator.   If  the  type  is  1  the  code  is  as 

follows 
1-X  coordinate  of  origin  (in  crystal  axes  system) 
2-Y  coordinate  of  origin 
3-Z  coordinate  of  origin 
4-to  angle  for  transformation,  to  standard,  orthonormal 

coordinate  system 
5~x  angle  for  transformation 
6-<J>  angle  for  transformation 

7-X  coordinate  of  an  atom  in  the  special  system 
8-Y  coordinate  of  an  atom 
9-Z  coordinate  of  an  atom 
10-T, 


11-T 
12 -T 
13-T 
14-T 
15-T 
16-L 
17-L 
18-L 
19-L 
20-L 
21-L 
22-S 
23-S 
24-S 
25-S 
26-S 
27-S 
28-S 
29-S 
30-S 


11 
22 
33 
12 
13 
23 
11 
22 
33' 
12 
13 
23 
11 
22 
33 
12 
13 
23 
21 
31 
32 


Columns  15-24   initial  value  of  this  parameter. 


11.   Extinction  and  Temperature  Factor  Type  Card: 
or  NEW(6)  i   0. 


Needed  only  if  NEW(5) 


Column  1  contains  the  extinction  factor  code: 

0  -  Isotropic,  do  not  change 

1  -  Anisotropic  Type  I  (mosaic  spread  dominates) 
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2  -  Anisotropic  Type  II  (domain  size  dominates) 

3  -  Now  Isotropic,  change  to  Type  I 

4  -  Now  Isotropic,  change  to  Type  II 

5  -  Now  Type  I,  change  to  Isotropic 

6  -  Now  Type  II,  change  to  Isotropic 

7  -  Now  Type  I,  change  to  Type  II 

8  -  Now  Type  II,  change  to  Type  I 

The  temperature  factor  type  code  is  punched  with  one  column  per  atom  be- 
ginning in  Col.  2  with  the  following  code: 

,0  -  Now  isotropic,  do  not  change 

1  -  Now  anisotropic,  do  not  change 

2  -  Now  isotropic,  change  to  anisotropic 

3  -  Now  anisotropic,  change  to  isotropic 

4  -  Now  3rd  cumulant,  do  not  change 

5  -  Now  anisotropic,  change  to  3rd  cumulant 

6  -  Now  3rd  cumulant ,  change  to  anisotropic 

7  -  Now  4th  cumulant ,  do  not  change 

8  -  Now  3rd  cumulant,  change  to  4th 

9  -  Now  4th  cumulant,  change  to  3rd 

12.  Parameter  Selection  Cards.   Needed  only  if  NEW(7)  4   0.   If  JXPAR  ±   0 
the  first  card  contains  refinement  selection  for  the  special  para- 
meters in  the  first  NPAR  columns.  A  zero  in  a  column  means  that  the 
corresponding  parameter  will  not  be  refined  and  a  one  means  that  it 
will."  After  this  cards  with  the  refinement  indicators  for  the 
conventional  parameters  punched  in  columns  1-80.   If  the  number  of 
atom  parameters  is  such  that  more  than  80  values  are  required,  the 
parameter  selection  information  should  be  continued  in  column  1  of 
the  next  card.  The  order  of  parameters  is: 

a)  occupancy  -  If  single  occupancy,  the  total  occupancy  is  refined. 
If  multiple  occupancy,  the  occupancy  of  species  A  is  refined. 

b)  x  coordinate 

c)  y_  coordinate 

d)  z_  coordinate 

e)  thermal  parameters,  one  for  isotropic  atoms,  six  for  anisotropic 
atoms.   The  thermal  coefficients  occur  in  the  same  order  as  on 
the  Atom  Parameter  Card  2 . 

f )  3rd  and  4th  cumulant  parameters  (if  needed) . 

13.  Parameter  Constraint  Cards (s)  :  Needed  only  if  NEW(7)  i   0. 
a)   Columns  Description 

1-3       Integer  specifying  the  dependent  atom  number. 

5-6       Dependent  parameter  number.  This  number  is  1  for  occup- 
ancy, 2-4  for  x  through  z  respectively,  5  for  B^  or 
Bn,  and  6,  7,  8,  9,  10  for  322>  833 »  &i2»  ^13  and 
3„_  respectively.   Values  of  11  to  20  treat  3rd 
cumulant  parameters  and  21-35  are  for  the  4th  cumulants, 
24 


13.  Parameter  Constraint  Cards(s)  (Cont.): 

a)  Columns  Description 

7-9       Atom  number  of  first  independent  parameter.   Special 

parameters  are  arbitrarily  indicated  as  atom  number  0. 

10-12      Parameter  number  of  first  independent  parameter.  Atom  0 
parameter  0  designates  an  initializing  constant  in  a 
linear  constraint  function. 

13-20      Coefficient  of  first  relationship. 

21-23      Atom  number  of  second  independent  parameter 

24-26      Second  independent 'parameter  number 

27-34      Second  coefficient 

35-37      Third  independent  atom  number 

38-40      Third  independent  parameter  number. 

41-48      Third  coefficient 

49-51      Fourth  independent  atom  number 

52-54      Fourth  independent  parameter  number 

55-62      Fourth  coefficient 

63-65      Fifth  independent  atom  number 

66-68      Fifth  independent  parameter  number 

69-76      Fifth  coefficient. 

If  the  independent  atom  number  is  zero,  designating  a  special  parameter, 
and  the  constraints  are  non  linear,  the  second  and  fourth  independent 
parameter  numbers  may  be  negative  to  denote  the  end  of  an  inclusive  list 
beginning  with  the  first  (third)  independent  parameter  number. 

b)  One  (1)  blank  card  or  a  card  with  zero  in  column  3. 

14.  Structure  Factor  Cards:  Needed  only  if  NEW(8)  4   0. 

a)  Format  Card:   This  card  contains  the  format  of  the  input  structure 
factor  data.  The  number  and  order  of  variables  on  a  card  (or  tape 
record)  depends  on  the  value  of  NEW(8) ,  as  explained  in  detail  below. 
All  records  must  contain  at  least  four  integer  fields,  for  the  indices 
and  an  observed/unobserved  indicator,  and  two  floating  point  fields 
for  F0bs  and  ap  or  other  weight  information.  They  may  contain 
other  information,  such  as  a  scale  factor  group  identifier,  the  factor 
6(6)  used  in  extinction  corrections,  and  the  setting  angles,  needed 
for  anisotropic  extinction  corrections. 

b)  Only  if  IWT=5  include  a  card  with  the  "ignorance  factor"  punched  as 
a  real  number  in  the  first  8  columns. 

c)  Only  if  NEW(8)=3  include  the  wavelength  in  A  in  the  first  8  columns. 
For  neutron  data  only  an  approximation  to  6(9)  will  be  computed  ac- 
cording to  the  formuja  $(0)  =  [ (2X2/(V2sin20) ]  x  108,  where  V  is  the 
unit  cell  volume  in  A3.   The  refined  quantity  then  will  be  the  product 
Tr*,  with  T  assumed  to  be  approximately  constant. 
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d)  Data  Cards  in  a  form  consistent  with  Format  Card: 
list  is  in  this  order. 


If  NEW(8)=1  the 


Variable  Description 

IH(1)      h  for  reflection 

IH(2)      k 

IH(3)       1 

FOBS        Observed  structure  factor 

SIGMA       Weight  information  consistent  with  IWT  on  Control  Card 

EXBETA     Factor  g(6)  for  extinction 

I  Number  of  scale  factor  for  this  reflection.  If  zero,  the 

scale  factor  is  the  same  as  the  preceding  reflection. 

Scale  factor  1  is  assumed  for  the  first  reflection. 
MREJ       Minimum  Observable  Parameter 

1  -  Intensity  above  minimum  observable. 

2  -  Intensity  below  minimum  observable. 
IEF       End-of-list  marker,  zero  on  all  data  cards 

TTH  20  These  angles  are  needed  only  if  anisotropic 

OMEG  co  extinction  is  to  be  calculated  and  must 

CHI  X  conform  to  the  diffractometer  conventions 

PHI  $  given  by  Busing  and  Levy  (1967,  p.  458). 

If  NEW(8)  =  2  the  input  list  is  XH(1),  IH(2),  IH(3),  MREJ ,1 , FOBS , EXBETA , 
SIGMA,  with  the  same  meanings  as  above.   The  end  of  the  data  must  be 
marked  by  an  end-of-file. 

If  NEW(8)  =  3  the  list  is  IH(1),  IH(2),  IH(3),  MREJ, FOBS, SIGMA. 
EXBETA  is  computed,  and  again  the  end  of  data  must  be  marked  by  an 
end-of-file. 

e)   End-of-list  card  -  A  card  with  a  non-zero  value  in  the  field  of 
variable  IEF  or  an  end-of-file  card  appropriate  to  the  system. 

15.  Run  Continuation  and  Parameter  Punch  Card:  Needed  for  all  cycles. 


Column 


Variable 


ICONT 


3-4 


JPUN 


ISC1 


Description 

-  No  more  cycles .  The  data  set  on  file  ISCI 
will  be  updated  if  IFILE  from  the  Control 
Card  is  zero 

-  Read  a  new  title  card  and  start  another  cycle 

-  Do  not  punch  any  atom  parameters 

-  Punch  atom  cards  in  form  required  by  'X-RAY67' 

-  Punch  atom  cards  in  form  required  by  this 
program 

3  -  Punch  atom  cards  in  both  forms 

If  this  field  is  non-zero  its  value  will  be  sub- 
stituted as  the  unit  number  for  the  binary 
master  file,  making  it  possible  to  try  new 
courses  of  refinement  without  destroying 
previous  stages . 
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BONDAN:   A  program  for  computing  bond  distances,  angles, 
and  thermal  ellipsoids,  with  errors 

This  program  uses  the  master  binary  file  created  by  RFINE  to 
compute  bond  distances  and  angles,  and  also  thermal  ellipsoids,  with 
associated  standard  deviations.   If  the  rigid  body  model  is  used, 
the  atomic  positions  of  the  rigid  body  will  be  corrected  for  lib- 
ration  before  calculation  of  distances  and  angles.   The  program  will 
also  accept  input  from  cards  by  means  of  the  program  BADTEA.   The 
complete  description  of  BADTEA  is  included  in  Appendix  B.   The 
details  of  data  input  are  as  follows: 

1.  Input  format  indicator  card 

Col.  1-2       1-input  from  cards  via  BADTEA. 

2-input  from  binary  file  created  by  RFINE. 

If  the  input  format  indicator  is  1,  see  appendix  B  for  further  input. 
If  the  input  format  indicator  is  2 ,  the  card  input  continues  as 
follows . 

2.  Title  Card:   Any  desired  hollerith  information  in  columns  1-72. 

3.  Search  Card:   Format (12 ,3A4) 

If  column  2  is  zero  or  blank,  a  search  of  the  data  tape  will  be  made 
to  find  a  data  set  with  a  label  the  same  as  columns  3-14.   The  error  cal- 
culation will  be  performed  on  that  data  set.   If  column  2  is  non-zero, 
the  first  data  set  on  the  tape  will  be  used  in  the  error  calculation. 

4.  Parameter  Card:   Format(2Il,3F8.0) 

Column     Variable  Description 

1  MBODAN      Bond  distance,  angle  selector 

0  -  None  computed 

1  -  Distances  only 

2  -  Distances  and  angles  computed 

2  MELVIB      Thermal  ellipsoid  selector 

0  -  No  ellipsoids  computed 

1  -  Thermal  ellipsoids  computed  for  anisotropic 

atoms 
3-10       DLIMIT(l)   Minimum  bond  distance  output 
11-18       DLIMIT(2)   Maximum  bond  distance  used  in  angles  output 
19-26       DLIMIT(3)   Maximum  bond  distance  output 

5.  Cell  variance-covariance  matrix:   Format (6E12 .6) 

This  matrix  is  the  matrix  output  by  the  lattice  constant  program  of 
C.  W.  Burnham. 
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Column  Description 

Card  1 

1-12  Variance  of  _a  (in  angstroms**2) 

13-24  Covariance  of  a_  with  b 

25-36  Covariance  of  a_  with  c 

37-48  Covariance  of  a_  with  alpha  (in  angstrom-radians) 

49-60  Covariance  of  a_  with  beta 

61-72  Covariance  of  a_  with  gamma 
Card  2 

1-12  Variance  of  b 

13-24  Covariance  of  b_  with  c_ 

25-36  Covariance  of  b_  with  alpha 

37-48  Covariance  of  b_  with  beta 

49-60  Covariance  of  b_  with  gamma 

61-72  Variance  of  c_ 
Card  3 

1-12  Covariance  of  £  with  alpha 

13-24  Covariance  of  c  with  beta 

25-36  Covariance  of  c   with  gamma 

37-48  Variance  of  alpha  (in  radians**2) 

49-60  Covariance  of  alpha  with  beta 

61-72  Covariance  of  alpha  with  gamma 
Card  4 

1-12  Variance  of  beta 

13-24  Covariance  of  beta  with  gamma 

25-36  Variance  of  gamma 

If  the  errors  in  the  cell  constants  have  a  small  effect  on  the 
distances,  the  last  4  cards  may  be  replaced  by  blanks. 
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FOURIER 

This  program  uses  the  binary  structure  factor  file  created  by  RFINE 
(ISC4-unit  12)  as  input  for  a  Fourier  synthesis.  The  map  is  written  on 
unit  11  for -input  to  the  plotting  program  CNTPLT  or  the  arbitrary  section 
program  ARBSECT.  The  details  of  data  input  are  as  follows: 

1.  Title  Card:  Any  desired  •alphanumeric  information  in  columns 
1-72.   This  information  heads  all  map  sections. 

2.  Parameter  Card:   Format  (611,  2X,  914,  3F8  0.4) 


Columns 


2. 


9-12 


Variable 


ITYPE 


ISECT 


INCENT 


KREJ 


IPLT 


LSPC 


IX 


13-16 

IXI 

17-20 

IXF 

21-24 

IY 

25-28 

IYI 

29-32 

IYF 

33-36 

IZ 

37-40 

IZI 

41-44 

IZF 

45-52 

EXXAG 

53-60 

F000 

Description 

Map  type  indicator 

0  -  Fourier  map 

1  -  Difference  Fourier 
Section  parameter 

0  -  x  across ,  v_  down  and  z_  sections 

1  -  x   " 


£ 
Z 

z 


X 

£ 

X 


£ 

x_ 
x_ 

£ 


Centric  indicator 

0  -  Centric 

1  -  Acentric 

0  -  Ignore  rejection  flag  and  include 

all  reflections. 

1  -  Do  not  include  any  reflections 

with  the  rejection  flag  set. 
If  non-zero,  map  sections  will  be 
written  to  unit  S2  for  plotting  or 
calculation  of  an  arbitrary  section. 
If  zero,  the  lines  of  the  map  will  be 
double  spaced.  A  non-zero  value  indi- 
cates single  spacing. 
Number  of  divisions  along  horizontal 
axis  (>0) 
Initial  point 
Final  point 

Number  of  division  along  vertical  axis 
(>0) 

Initial  line 
Final  line 

Number  of  divisions  along  section  axis 
(>0) 

Initial  section 
Final  section 
Map  exaggeration  desired 
Number  of  electrons  in  unit  cell  for 
Fourier . 
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2.  Parameter  Card  (cont.): 

Zero  for  difference  Fourier. 
61-68        TMULT       Multiplicity  of  a  general  reflection 

3.  Multiplicity  Card:   Format  (12F6.6 


Columns 

Description 

1-6 

Multiplicity  of  hOO  or  OkO 

7-12 

Multiplicity  of  001 

13-18 

Multiplicity  of  hkO 

19-24 

Multiplicity  of  hOl 

25-30 

Multiplicity  of  Okl 

31-36 

Multiplicity  of  hhO 

37-42 

Multiplicity  of  hOh 

43-48 

Multiplicity  of  Okk 

49-54 

Multiplicity  of  hhh 

55-60 

Multiplicity  of  hhl 

61-66 

Multiplicity  of  hkk  or  hkh 

67-72 

Multiplicity  of  hkl 

C.  Size  of  Program: 

This  program  uses  the  array  FBUF  in  common  block  ABLE  as  pooled 
storage.   The  size  of  this  array,  S,  is  governed  by  the  map  to  be  cal- 
culated and  may  be  estimated  using: 

S  =  3N  ,  +  2(h'    +  2)  (2k*    +1) 
ref       max        max 

+  2(h*    +  l)(Pv  +  1) 
max       x 

where  Nref  in  the  number  of  reflections  included,  h'max  is  the  maximum 
value  of  the  index  associated  with  the  horizontal  axis,  k'   x  is  the 
maximum  value  of  the  index  associated  with  the  axis  plotted  down  the  page, 
and  P  is  the  number  of  points  calculated  in  the  horizontal  direction. 
If  the  Fourier  sections  are  to  be  saved  for  later  plotting,  S  should 
be  increased  by  P  P  .   If  the  computer  operating  system  supports  dynamic 
core  allocation,  the  program  may  be  converted  to  automatically  request 
the  proper  memory  size. 

D.  Output: 

Each  map  point  is  converted  to  a  four-digit  integer  before  print- 
ing with  zero  values  replaced  by  blanks.   Each  point  will  occupy  an  area 
4  columns  by  1  or  2  rows  depending  upon  the  spacing  selected.   For  most 
printers,  the  area  is  0.4  x  1/6  inches  for  single  spacing  or  0.4  x  1/3 
inches  for  double  spacing  of  the  lines.   If  the  map  requires  more  than 
30  points,  the  section  will  be  printed  in  strips  with  30  points  per  strip. 

30 


E.  Additional  Problems: 

If  additional  maps  are  to  be  run  using  the  same  structure  factor 
list,  a  complete  set  of  cards  beginning  with  the  Title  Card  should  be 
included.   An  end-of-file  in  the  card  stream  terminates  the  program. 
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LISTFC 

This  program  is  an  adaption  of  the  XRAY67  program  LIST  FC  (Stewart 
et  al.,  1967).   It  has  been  altered  to  stand  alone  and  read  the  structure 
factor  file  (ISC4  -  unit  12)  created  by  RFINE.   The  output  of  the  program 
has  been  changed  so  that  the  "phase  in  millicycles"  has  been  replaced  by 
sigma  (Fobs).   Input  is  as  follows: 

Card  1. 

Col.  1-6  compound  identification;  Col.  7-78  any  title. 
Card  2. 

Col.  14        (l)/(2)/(3)  for  h  index  varies  (most) /(next  most)/ 
(least)  rapidly 

15  Same  for  k  index 

16  Same  for  I   index 

(Note  that  the  sum  and  product  of  Cols.  14,  15,  and  16 
must  equal  6.) 

17-20      Number  of  lines  per  list  page. 

21        (blank) /(A)  for  (DO) /(DO  NOT)  internally  change  the 
value  of  the  number  of  lines  in  order  to  make  the 
bottom  of  the  last  page  as  even  as  possible. 

22-24      Number  of  list  columns  per  page.  Note  that  the  pro- 
duct of  the  number  of  lines  and  the  number  of 
columns  cannot  exceed  2000 .   That  is  no  more 
than  3000  items  for  the  number  of  reflections  and 
headings  (including  spaces)  per  page. 

26        Number  of  blank  print  columns  before  the  LISTFC 
column. 

28        Number  of  print  columns  for  most  rapidly  changing 
index. 

30        Number  of  print  columns  for  F0^s 

32        Number  of  print  columns  for  Fca^c 

34        Number  of  print  columns  for  ap 

36        Special  flag  for  unobserved  reflections  (blank=*) 

38        Special  flag  for  severely  extinct  reflections 
(blank=E) 

40        Special  flag  for  special  reflection  (JC0DE=4) 

The  next  eight  fields  are  (blank) /(l)  for  (do  not) /do) 

42        Print  symbol  for  unobserved  or  severely  extinct 
reflections. 
This  symbol  switch  adds  one  print  column 

44        For  centric  structures,  attach  sign  of  A  to  F  . 

46        Double  space  the  lines 

48        Restore  each  LISTFC  page  to  the  top  of  the  printer 

page 
50        Print  title  at  the  top  of  each  printer  page 
52        Punch  a  set  of  "FCARD"  cards 

56        Print  a  minus  sign  on  F0bs  for  unobserved  reflections 
58        Write  a  separate  copy  of  FC  list  on  unit  7 
62        Number  of  times  to  try  to  get  headings  at  the  top  of 

all  columns 
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64         (l)/(2)/(3)  for  headings  separated  by  (no  blank 
lines) /(l  blank  line  above) /(l  above  and  1 
below) (blank=3) 
65-72      factor  Fs  are  to  be  multiplied  by  (blank=10) 
73-80      factor  cts  are  to  be  multiplied  by  (blank=10) 

Note  that  all  F  and  a  values  are  printed  as  integers, 
so  that  scale  factor  must  be  sufficient  to  preserve  a 
sufficient  number  of  significant  figures . 


33 


CNTPLT 

CNTPLT  is  a  contour  plotting  program  for  Fourier  maps  and  other 
similar  functions.   The  routines  contained  in  this  program  have  been 
derived  from  multiple  sources  with  the  original  version  of  subroutine 
ICONT  written  by  S.  H.  Zisk  and  N.  M.  Brenner  of  M.I.T.  Lincoln  Lab. 
Subroutine  INTPLT  uses  the  rational  polynominal  spline  interpolation 
of  D.  V.  Ahujo,  I.B.M.  Systems  Journal,  p.  208-217,  (1968).   The  other 
routines  have  been  modified  from  code  prepared  by  G.  Ford,  Purdue  Univer- 
sity. 

The  present  version  of  these  routines  have  been  written  to  be  as' 
machine  independent  as  possible.   It  is  believed  that  only  the  call  to 
the  plotter  function  PLOTS  (Sequence  No.  A5300)  and  the  two  functions 
MARK  and  UNMARK  called  by  ICONT  are  machine  specific.   This  is  the  main 
contribution  of  the  present  author  (LWF) . 

Method 

The  main  contouring  subroutine,  ICONT,  takes  the  function  values 
at  the  grid  points  of  the  map  and  computes  the  coordinates  of  the  inter- 
section of  a  specified  contour  line  with  the  lines  of  the  grid.  Linear 
interpolation  is  used  in  this  process.  The  second  contouring  subroutine, 
CPLOT,  coordinates  the  drawing  of  the  contours  and  produces  any  desired 
contour  labels.   The  third  contouring  routine,  INTPLT,  produces  a  smooth 
curve  through  the  intersections  located  by  ICONT.  Rational  cubic  poly- 
nomials are  used  to  fit  a  spline-like  curve  through  neighboring  points . 
Two  such  curves  meet  smoothly  because  the  first  and  second  derivatives 
are  single  valued  along  the  composite  curve.  The  plotter  line  is  drawn 
as  a  series  of  straight  line  segments  with  the  number  of  such  segments 
specified  by  the  user. 

The  main  program  includes  various  options  which  allow  labeling  of 
plot  sections  and/or  grid  points,  drawing  of  the  bounding  parallelogram, 
and  other  similar  functions. 

Input 

The  values  of  the  function  to  be  plotted  are  read  from  a  binary 
input  file  (unit  10)  with  each  section  occupying  two  records.   The 
first  contains  the  variables  MA,  MB,  JZ,  MZ  where  MA  and  MB  are  the 
number  of  grid  points  along  x  and  y  respectively.   The  second  record 
contains  the  data  for  the  map  section  at  JZ/MZ.   These  data  are  assumed 
to  be  in  real  form  and  are  packed  in  an  array  which  is  MAxMB.   The 
output  file  produced  by  program  FOURIER  is  in  this  form. 

Input 

1.  Title  Card:   Format  (20A4) 

2.  Control  Card:   Format  (F5.0,  15 ,3F5. 0,615, F5. 0,215) 
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Columns 


Variable 


Description 


1-5 

6-10 
11-15 


16-20 


21-25 


TLEHGT 

NLABL 
A 


ALP 


Plotting  height  of  title  in  inches.   If 

<0,  no  title  will  be  plotted 

Number  of  vertex  labels  to  be  plotted 

Width  of  plot  (in- inches)  parallel  to 

x-axis  of  Fourier.   If  A>0,  an  11' 

plotter  is  assumed.   If  a  30"  machine 

is  available,  A  should  be  less  than  zero. 

Length  of  plot  (in  inches)  parallel  to 

v_-axis  of  Fourier  which  is  parallel  to 

plotter  x-axis. 

Angle  between  Fourier  axes  (in  degrees) 


26-30 

NA 

31-35 

NB 

36-40 

NP 

41-45 

NCONT 

46-50 

NLINE 

51-55 

NCROSS 

56-60 

SCROSS 

61-65 

NSECT 

66-70 

NZ 

Number  of  grid  points  parallel  to  Fourier 

x_-axis 

Number  of  grid  points  parallel  to  Fourier 

v_-axis 

Number  of  interpolation  points  per  spline. 

The  default  value  is  4 . 

Number  of  sets  of  contour  levels  (NCONT 

£10) 

Number  of  line  vertices  specified 

(NLINE<20) 

Number  of  crosses  to  be  placed  on  plot 

(NCROSS<20) 

Size  of  cross  (in. inches) 

Number  of  sets  of  sections  specified 

Grid  spacing  of  section  levels. 


3.   Contour  Card(s) :   Format  (3F5.0,I5) 
There  should  be  NCONT  of  these  cards, 


Columns 

1-5 

6-10 

11-15 

16-20 


Variable 

SBEG 
CSTEP 
SLIM 
NSIG 


Description 

Height  of  initial  contour 

Step  between  contours 

Height  of  final  contour 

Number  of  digits  after  decimal  point  in 

contour  label.   If  NSIG<0,  contours  will 

not  be  labelled. 


4.  Line  Vertex  Card (s) :   Format  (16F5.0) 

The  positions  of  the  line  vertices  in  map  grid  coordinates  should 
be  punched.  A  line  will  be  drawn  from  point  i_  to  point  i+1  if  x^+,HD. 
No  line  will  be  drawn  if  xi+-i<0.   Omit  these  cards  if  NLINE<0. 
Note:   The  first  grid  point  is  (1,1),  not  (0,0).   i.e.  add  (1,1)  to 
scaled  coordinates. 

5.  Cross  Card(s):   Format  (16F5.0) 
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The  map  grid  coordinates  of  the  crosses  should  be  punched  sequentially. 
Omit  these  cards  if  NCROSS<0. 

6.   Section  Selection  Cards:   Format  (315) 

Columns        Variable        Description 


1-5 

6-10 

11-15 


INITZ 

ISTEPZ 

IFINZ 


Sections  from  INITZ 

to  IFINZ  in  steps 

of  ISTEPZ  will  selected  and  contoured. 


There  should  be  NSECT  of  these  cards. 

7.   Vertex  Label  Card (s) :   Format  (5F5.0,A6) 


Columns 

1-5 
6-10 
11-15 
16-20 
21-25 
26-31 


Variable 

WX 

WY 

DX 

DY 

HGT 

LABL 


Description 

x-grid  coordinate 
y-grid  coordinate 
Plotter  x-axis  offset 
Plotter  y-axis  offset 
Label  height  in  inches 
Alphanumeric  label 


There  should  be  NLABL  of  these  cards . 

Tape  output  for  Calcomp  plotter  is  on  Unit  12, 
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ARBSECT 


This  program  calculates  an  arbitrary  section  of  a  Fourier  map  from 
a  complete  set  of  sections  calculated  in  regular  intervals  parallel  to 
the  crystallographic  axes.   The  input  sections  may  be  computed  with  the 
regular  Fourier  program  and  are  assumed  to  be  in  file  No.  11.   The  plane 
of  the  output  section  is  specified  by  the  coordinates  of  the  plane  to 
be  calculated. 

Program  ARBSECT  has  several  options  available.   If  only  a  fraction 
of  the  unit  cell  is  used  as  input,  it 'will  calculate  the  electron  density 
in  the  accessible  portion  only  with  the  remainder  replaced  by  asterisks 
in  the  printout.  However,  if  the  map  for  a  complete  unit-translation  has 
been  computed,  the  program  will  use  the  translation  symmetry  to  calculate 
a  complete  section.   No  other  symmetry  operations  are  used.   Other  options 
include  (a)  describing  the  points  used  to  define  the  plane  in  terms  of 
map  or  fractional  coordinates  as  desired,  (b)  saving  the  resulting 
section  in  file  12  for  later  plotting,  (c)  computing  the  entire  section 
or  selecting  a  portion  of  it,  and  (d)  calculating  and  printing  the  geometry 
including  axial  lengths  and  interaxial  angle,  of  the  output  section. 

Method 

The  input  coordinates  of  the  three  points  defining  the  plane  are 
initially  converted  to  map  coordinates  if  necessary.  The  coefficients 
of  the  plane  in  the  equation  a,x  +  a,  y  +  a,z  +  a,  =  0  may  be  found  from 
the  determinant  equation 

x    y    z    1 

xl  yl  zl  X 
x2  y2  Z2  1 
x3   y3   Z3   1 

where  x^,  y^,  and  z^  are  the  coordinates  of  the  i'th  point.  This  plane 
will  have  a  line  of  intersection  with  each  section  of  the  original  map. 
For  the  values  of  y_  and  z_  corresponding  to  each  line  and  section  of  the 
output  map,  the  x-coordinate  of  this  point  in  the  plane  is  determined. 
A  cubic  function  of  the  form  ax3  +  bx2  +  ex  +  d  is  fitted  to  the  four 
adjacent  grid  points  on  the  line,  and  the  value  of  this  function  at  the 
intersection  point  is  used  for  the  output  at  that  y  and  z.   This  process 
is  repeated  until  all  output  points  have  been  determined. 

Input 

The  input  map  sections  are  read  from  file  11  and  the  resulting  map 
is  output  to  file  12.   These  assignments  may  be  changed  by  redefining 
variables  ISC  and  ISC1  respectively. 
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1.  Title  Card  -  Format  (18A4) 

Columns  1-72  of  this  card  may  contain  any  desired  information. 

2.  Control  Card  -  Format  (813) 


columns 


1-3 


Variable 


10 


6 
9 

12 


IZO 
ISPC 

IMRK 


13-15 

IYI 

16-18 

IYF 

19-21 

IZI 

22-24 

IZF 

Description 

If  IC#0,  the  output  section  will  be 
written  on  Unit  12. 
Section  number  of  output 
If  ISPC  =  0  the  output  map  will  be 
single  spaced;  double  spaced  if  ISPC=1. 
If  IMRK=0,  the  y-  and  z-  limits  of  the 
input  sections  will  be  used  as  the  x-, 
and  y-limits,  respectively,  of  the  output 
section.   If  IMRK  ^0,  "these  values  are 
read  from  this  card. 

Initial  value  of  x  in  transformed  section 
Final  value  of  x  in  transformed  section 
Initial  value  of  y  in  transformed  section 
Final  value  of  y  in  transformed  section 


3.   Cell  Card  -  Format  (6F8.0) 


Columns     Description 

o 

1-8  Length  (in  A)  of  cell  axis  across  page  in  original  section 

9-16  Length  of  cell  axis  down  page  in  original  section 

17-24  Length  of  cell  axis  out  of  page  in  original  section 

25-32  Angle  (in  degrees)  between  axis  1  and  2 

33-40  Angle  between  axes  1  and  3 

41-48  Angle  between  axes  2  and  3 

The  above  length  and  axes  are  required  only  if  the  program  is  to 
compute  and  print  the  axial  length  and  interaxial  angle  of  the  transformed 
section.   If  not,  this  card  may  be  blank. 

4.   Plane  Card  -  Format  (9F8.0,  12) 


Columns 

Description 

1-8 

x-  coordinate 

of   first  point 

9-16 

y- 

ii        it            ii 

17-24 
25-32 

z-            " 
x-            " 

ii        it            ii 
,,  second      „ 

33-40 

y- 

ii        ii            ii 

41-48 
49-56 

z-             " 
x-             " 

ii        ii            it 
„  third        „ 

57-64 

y-         " 

ii        ii            ti 
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4.   Plane  Card  -  (Cont.) 

Columns     Description 

65-72      z-  coordinate  of  third  point 

74        Coordinate  type  indicator.   If  zero,  the  coordinates  are 

given  in  grid  spacings.   If  non-zero,  the  points  are 

given  in  fractional  coordinates. 

Additional  arbitrary  slices  may  be  computed  by  supplying  additional 
data  decks  beginning  with  the  Title  Card. 

Form  of  Input  Sections 

The  input  Fourier  sections  are  read  from  Unit  11  with  unformatted 
reads.   There  must  be  two  records  per  section  as  follows: 


Record  1: 

Word 

Variable 

1 

JXA 

2 

JYA 

3 

KZC 

4 

IZ 

5 

IX 

6 

IXI 

7 

IXF 

8 

IY 

9 

JYI 

10 

JYF 

11 

JZI 

12 

JZF 

Record  2 : 

Description 

Number  of  grid  points  calculated  parallel  to  x 

Number  of  grid  points  calculated  parallel  to  y 

Section  number 

Number  of  sections  per  unit  cell 

Number  of  intervals  per  unit  cell  parallel  to  x 

Initial  grid  points  parallel  to  x 

Final  grid  point  parallel  to  x 

Number  of  intervals  per  unit  cell  parallel  to  y 

Initial  grid  point  parallel  to  y 

Final  grid  point  parallel  to  y 

Initial  grid  point  parallel  to  z 

Final  grid  point  parallel  to  z 


This  record  consists  of  the  JXA*JYA  values  of  the  map  in  this 
section.   These  real  values  must  be  packed  in  an  array  with  the  x  co- 
ordinate as  the  first  subscript  and  y  as  the  second  subscript. 
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Appendix  A 
Overlay  scheme  for  RFINE 

When  program  RFINE  is  run  under  operating  systems  which  charge 
according  to  the  product  of  core  allocation  and  CPU  time,  substantial 
savings  can  be  achieved  by  using  an  overlay  structure  which  causes  only 
that  portion  of  the  program  which  is  actually  needed  to  be  resident  in 
core.   Furthermore ,  there  can  be  large  savings  if  the  core  allocated 
for  the  dependency  list  and  the  normal  equations  matrix  can  be  adjusted 
as  appropriate  to  a  particular  problem.   On  UNIVAC  machines  operating 
under  EXEC  8  this  latter  function  is  accomplished  dynamically  at  run 
time  by  means  of  the  subroutine  CORSZ,  which  makes  use  of  executive 
requests  ER  MCOREg  and  ER  LCOREg .  On  other  machines  which  don't  have 
a  similar  function  capability,  this  can  be  achieved  by  a  run  time 
recompilation  of  the  short  main  program,  with  the  dimension  of  the 
array  A  in  common  block  LSMAT  appropriately  adjusted. 

The  program  has  been  run  successfully  with  a  structure  consisting  of 
a  main  stem  and  three  overlay  segments  as  follows: 

Main  stem:  main  program,  subroutines  RFINE  and  CORSZ. 
Overlay  1:   Subroutines  RESET, DERIVS , RDEGN ,  and  PAREXT  plus 

auxilliary  functions  called  by  them,  including  DIJKLM.DETERM, 
EL4,PIJKLM,BIJK,CIJKL,AIJ,VXV,S0LVE,EIJ,FIJ,IPAK,PDEGN,D3KDS, 
D2KDS,P3KDL,P2KDL,D2KDT,TKTLS,TFM,SKTLS,CORECT,MULMM,TRNSPS, 
EULER,  FILLIN,SPVAL,ERRCAL. 
Overlay  2:   Subroutines  INPUT,SFAC,RCALC,MATRIX,MODIFY,  and  SYMINV. 
Overlay  3:   Subroutines  TRNFRM , B0DAN1 ,  and  ELVIB1. 
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Appendix  B 
BADTEA-Alternate  input  for  bond  distances  a-nd  angles,  and  thermal  ellipsoids 

A.  IDENTIFICATION 

TITLE:   Crystallographic  error  analysis 

NAME:   U  M  BADTEA 

PROGRAMMER:   Larry  W.  Finger 

DATE:  December  27,  1965  (revised  November  1968) 

B.  PURPOSE 

U  M  BADTEA  (University  of  Minnesota  program  for  computing  Bond 
angles  and  Distances,  and  Thermal  Ellipsoids  with  error  Analysis) 
will  calculate  bond  distances,  bond  angles,  thermal  ellipsoids 
and  the  errors  associated  with  each  computed  value.   The  program 
was  converted  to  Fortran  IV  with  the  assistance  of  C.  T.  Prewitt. 

C.  USAGE 

1.  Input  deck  order: 

Title  card 

Cell  card 

Function  card 

Symmetry  card(s) 
*Standard  Deviations  of  cell  parameters 
Correlation  matrix  of  cell  parameters 

Atom  cards 
*Temperature  Factdr  Standard  Deviations 
*Temperature  Factor  Correlation  matrix 
*Atomic  Position  Standard  Deviations 
*Atomic  Position  Correlation  Matrix 

*These  data  may  not  be  needed.   The  next  section  describes  in 
detail  the  conditions  requiring  these  input  data.  All  input 
from  Atom  Cards  to  the  end  of  the  deck  may  come  from  a  mod- 
ified ORFFF  tape. 

2,  Description  of  input  deck: 

a.   Title  Card  FORMAT  (18A4) 

Columns  Description 

1-72  Identifier  for  problem.   This  should 

be  any  title  which  identifies  the 
compound . 
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b.   Cell  Card 


FORMAT  (6F10.5) 


Columns 

1-10 
11  -  20 

21  -  30 
31  -  40 
41  -  50 
51  -  60 


Description 

o 

a_  axis  length,  in  A 
b_  axis  length, 
c_  axis  length* 
a  in  degrees 


Function  Card 


Columns 

1-4 

5-8 


9-12 


13  -  16 


17  -  20 


24 


31  -  40 
41  -  50 
51  -  60 


d.   Symmetry  Cards 


Columns 

1  - 

3 

5  - 

6 

8  - 

9 

11  - 

13 

15  - 

16 

18  - 

19 

21  - 

23 

25  - 

26 

28  - 

29 

FORMAT  (514,  3X.A1.6X,  3F10.5) 

Description 

Number  of  atoms  in  asymmetric  unit  (<60) 
Centric  Indicator 

0  -  Structure  is  centric 

1  -  Structure  is  non-centric 
Bond  Distance,  Angle  Parameter 

0  -  No  angles  or  distances  calculated 

1  -  Distances  only 

2  -  Angles  and  Distances 
Ellipsoid  Indicator 

0  -  None  calculated 

1  -  Ellipsoids  calculated  for  aniso- 

tropic atoms 
Error  Indicator 

0  -  Errors  not  calculated 

1  -  Errors  calculated . 

2  -  Errors  calculated  with  structure 

information  read  from  modified 
ORFFE  tape. 
Lattice  Type  Code  -  This  column  should 
be  punched  with  a  P,  A,  B,  C,  F,  I,  or 
R  depending  on  the  lattice  type. 
Minimum  bond  distance  to  be  printed 
Maximum  distance  used  in  angle  calculation, 
Maximum  bond  distance  to  be  printed. 

FORMAT  (3 (A3, 2 (IX, A2), IX), 14) 

Description 

Translational  part  of  x_' 
First  positional  part  of  jc' 
Second  positional  part  of  x_' 
Translational  part  of  y_' 
First  positional  part  of  y_' 
Second  positional  part  of  y_' 
Translational  part  of  z_' 
First  positional  part  of  z_? 
Second  positional  part  of  a* 
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Columns 
34 


Description 

End  of  Symmetry  Deck  Indicator 

Blank  -  more  symmetry  cards  follow 
1  -  this  is  the  last  card 


Note:   The  translational  parts  of  the  symmetry  operators  are  of  the 
form:  bbb  1/2  1/3  2/3  1/4  3/4  1/6  or  5/6  where  b  signifies  a  blank 
The  positional  parts  are  of  the  form:  bw  +  w  -  w  where  w  may  be 
X,  Y  or  Z.   Symmetric  positions  which  are  related  by  a  center  of 
inversion  or  a  multiple  lattice  point  must  hot  be  included  since 
these  positions  are  internally  generated.   However,  position  x_  y_  z_ 
must  be  included  and  must  be  first. 

e.   Standard  Deviations  of  Cell  Parameters   FORMAT  (6F10.8) 

This  card  and  the  cards  in  section  f .  below  are  needed 
only  if  errors  are  to  be  calculated. 


Columns 

1-10 
11  -  20 
21  -  30 
31  -  40 
41  -  50 
51  -  60 


Description 

Standard  Deviation  of  a_ 
Standard  Deviation  of  b_ 
Standard  Deviation  of  c_ 
Standard  Deviation  of  a  in  degrees 
Standard  Deviation  of  3 
Standard  Deviation  of  y 


f.   Cell  Parameter  Correlation  Matrix   FORMAT  (6F8.5) 


Columns 

Corre 

dat  ion  C 

Card  1: 

1-8 

a 

with 

a_ 

(1.0) 

Card  2: 

1-8 

a 

with 

b 

9-16 

b 

with 

b 

(1,0) 

Card  3: 

1-8 

a 

with 

c_ 

9-16 

b 

with 

c_ 

17  -  24 

c_ 

with 

c^ 

(1.0) 

Card  4: 

1-8 

a 

with 

a 

9-16 

b 

with 

a 

17  -  24 

c_ 

with 

a 

25  -  32 

a 

with 

a 

(1.0) 

Card  5: 

1-8 

a 

with 

6 

9-16 

b 

with 

3 

17  -  24 

c 

with 

3 

25  -  32 

a 

with 

3 

33  -  40 

3  with 

3 

(1.0) 
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Cell  Parameter  Correlation  Matrix  (Cont) . 


Columns 

Correlatioi 

i  Coefficiei 

xt   of 

Card  6: 

1-8 

a   with  y 

9-16 

b_  with  Y 

17  -  24 

c_  with  Y 

25  -  32 

a  with  Y 

33  -  40 

3  with  Y 

41  -  48 

Y  with  Y 

(1.0) 

Atom  Cards 

FORMAT  (2A3, 

7X,  3F7.5, 

I2/6F8.5) 

There  should  be  a  set  of  2  atom  cards  for  each  atom 
in  the  asymmetric  unit. 

Position  Card: 

Columns  Description 

1-6  Identification  for  atom 

14-20         x_  in  fractional  coordinates 
21-27         y_  in  fractional  coordinates 
28-34         z_  in  fractional  coordinates 
36  Temperature  factor  type 

0  -  isotropic 

1  -  anisotropic 

Temperature  factor  card : 

1-8  6  or  B 

9-16  B„ 

17  -  24  3„ 

25  -  32  B?? 

33  -  40  B7, 

41  -  48  0g 

The  anisotropic  coefficients  should  be  those  for  which  the 
temperature  factor  expression  is 

exp  [-(h2$u  +  k2322  +  12333  +  2hk312  +  2hl313  +  2kl323)]. 

h.   Temperature  Factor  Standard  Deviations   FORMAT  (6F10.8) 

These  cards  are  needed  only  for  anisotropic  atoms  when 
ellipsoids  and  errors  are  to  be  calculated. 

Columns  Description 

1-10  Standard  Deviation  of  B,-, 

11  -  20  Standard  Deviation  of  B22 

21  -  30  Standard  Deviation  of  B„ 
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Temperature  Factor  Standard  Deviations  (Cont) 

Columns  Description 

31  -  40  Standard  Deviation  of  g_2 

41  -  50  Standard  Deviation  of  g-„ 

51  -  60  Standard  Deviation  of  Rl" 

i.   Temperature  Factor  Correlation  Matrix   FORMAT  (6  F8.5) 

These  cards  are  needed  only  for  anisotropic  atoms  when 
ellipsoids  and  errors  are  to  be  calculated. 

Columns  Correlation  Coefficient  of 

Card  1: 
1-8 
Card  2 : 
1-8 
9-16 
Card  3: 
1-8 
9-16 
17  -  24 
Card  4: 
1-8 
9-16 
17  -  24 
25  -  32 
Card  5: 
1-8 
9-16 
17  -  24 
25  -  32 
33  -  40 
Card  6: 
1-8 
9-16 
17  -  24 
25  -  32 
33  -  40 
41  -  48 

Each  anisotropic  atom  has  a  Standard  Deviations  card  and  a 
Correlation  Matrix. 

j.   Standard  Deviations  of  Atomic  Positions   FORMAT  (3F10.8) 

These  cards  are  needed  only  if  bond  distances  or  angles 
are  to  be  calculated  with  errors . 
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;..  with 

hi 

(1.0) 

!- 1  with 
L-  with 

R22 

322 

(1.0) 

i- .   with 

?33 

\\   with 
;33  With 

B33 
333 

(1.0) 

I-,  with 

ol2 

;  with 

2" 

oo  with 

P^  n 

12  with 

ft12 
312 

(1.0) 

, ,  with 

3i  0 

11 

ol3 

22  with 

ol3 

oo  with 

^3 

.a  with 

37o 

*3  with 

13 

313 

(1.0) 

with 

Boo 

o23 

2  with 

??3 

with  g" 
it  with  3„ 

zt   with  " 

Po  0 

23  With 

R23 
323 

(1.0) 

Standard  Deviations  of  Atomic  Positions  (Cont) 

Columns  Description 

1-10  Standard  deviation  of  x 

11  -  20  Standard  deviation  of  y_ 

21  -  30  Standard  deviation  of  z_ 

There  should  be  one  card  of  this  type  for  each  atom  in  the 
input  list. 

k.   Atomic  Positions  Correlation  Matrix    FORMAT  (10F8.5) 

These  cards  are  needed  only  if  bond  distances  or  angles  are 
to  be  calculated  with  errors.   The  lower  triangle  of  this 
matrix  is  read  by  rows  with > each  row  starting  on  a  new 
card.  When  a  row  contains  more  than  ten  elements,  it  should 
be  continued  on  the  next  'card. 

Correlation  Coefficient  of 

x_  of  atom  1  with  x_  of  atom  1    (1.0) 

x  of  atom  1  with  y_  of  atom  1 

y_  of  atom  1  with  y_  of  atom  1    (1.0) 

x_  of  atom  1  with  z^  of  atom  1 
y_  of  atom  1  with  x_  of  atom  2 
z_  of  atom  1  with  z_  of  atom  1    (1.0) 

x  of  atom  1  with  x  of  atom  2 

y_  of  atom  1  with  x.  of  atom  2 

x   of  atom  1  with  x  of  atom  2 

x_  of  atom  2  with  x  of  atom  2    (1.0) 

This  scheme  should  be  continued  until  all  elements  of  the 
matrix  have  been  punched. 


Columns 

Card 

1: 

1 

-  8 

Card 

2: 

1  - 

-  8 

9  ■ 

-  16 

Card 

3: 

1  ■ 

-  8 

9  - 

-  16 

17  ■ 

-  24 

Card 

4: 

1  - 

-  8 

9  - 

-  16 

17  - 

-  24 

25  - 

-  32 

D.   METHOD 


Bond  distances : 

If  X  is  a  column  matrix  of  fractional  atomic  coordinates  for  an 
atom  in  the  asymmetric  unit  and  Y  is  the  column  matrix  of  the 
coordinates  of  a  symmetry  related  atom,  then  Y  and  X  satisfy 
the  following  equation: 

Y  =  RX  +  T  (1) 

where  R  is  a  rotation-like  operation  in  that  it  involves  axis 
interchanges  and  T  is  a  translation  operation  which  includes 
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glide  and  screw  components  as  well  as  origin  shifts.  Defining 
Z  as  a  column  matrix  whose  coefficients  are  the  differences  in 
fractional  coordinates  between  an  atom  at  X  and  an  atom  at  Y, 
then 

Z  =  X  -  Y  (2) 

If  G  is  defined  as  the  real  space  metric  tensor  which  has  elements 

gij  =  ^  .  a.  (3) 

where  the  a^ s  are  the  real  space  axis  vectors,  then  G  is  the 
following  matrix: 

a     ab  cos  y 

G  =(  ab  cos  y    °  be  cos  a  )  (4) 

,ac  cos  3  be  cos  a 

T 
If  Z  denotes  the  transpose  of  Z  which  involves  an  interchange  of 

the  rows  and  columns  of  Z ,  it  may  be  shown  that  the  square  of  the 
bond  distance,  s^,  between  the  atoms  is  determined  by  the  fol- 
lowing matrix  equation: 

2     T 
s  =  Z  GZ  (5) 

In  the  program,  all  atoms  within  a  spherical  shell  centered  on 
each  atom  in  the  asymmetric  unit  are  located  and  the  distances 
are  calculated.   The  size  of  the  spherical  shell  is  determined 
by  the  minimum  and  maximum  distances  to  be  considered  which  are 
input  parameters.   The  search  for  the  contents  of  the  shell  is 
made  efficient  by  enclosing  the  central  atom  with  a  unit  cell 
shaped  parallelopiped  tangent  to  the  sphere  with  all  atomic 
locations  outside  this  parallelopiped  excluded.   The  position 
of  those  atoms  which  lie  within  the  spherical  shell  are  stored 
in  a  table.   This  table  also  includes  the  atom  number  and  the 
number  of  the  R  matrix  used  to  generate  this  position  which 
is  needed  to  properly  compute  the  partial  derivatives  of  the 
bond  distance  with  respect  to  the  positional  coordinates.   These 
two  numbers  are  stored  in  packed  form  in  the  same  storage 
location.   Each  time  a  new  atom  is  generated,  this  table  is 
searched  to  determine  whether  this  atom  has  already  been  included. 
This  search  assures  that  atoms  in  special  positions  are  properly 
treated  although  they  may  be  generated  in  several  different  ways . 

2 .   Bond  angles : 

All  possible  bond  angles  about  each  atom  in  the  asymmetric  unit 
are  computed  using  the  table  of  positions  prepared  during  the  bond 
distance  calculation.   Each  pair  of  atoms  from  this  table  is 
checked  to  determine  whether  their  distances  from  the  central  atom 
are  within  the  limits,  and  if  so,  the  bond  angle  is  computed.   If 
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X  denotes  the  column  matrix  of  the  central  atom  with  Y-,    and  Y2 
denoting  the  column  matrices  of  the  peripheral  atoms,  the  bond 
angle  0  for  the  angle  Y^  -  X  -  Y2  can  be  shown  to  satisfy  the 
following  relationship : 

(Y  -X)T  G(Y  -X) 

cos  0  =  ^ 5 (6) 

blX  b2X 

Where  G  is  the  real  space  metric  tensor  (see  equation  (4)) 

S.   is  the  bond  distance  from  the  atom  of  X  to  the 
atom  at  Y  and 

S„  is  the  bond  distance  from  the  atom  at  X  to  the 
atom  at  Y„ . 

3.   Ellipsoids  of  Vibration: 

The  conversion  of  anisotropic  thermal  coefficients  is  a  problem 
of  finding  the  characteristic  values  and  characteristic  vectors 
of  a  particular  matrix.   Note  that  characteristic  value,  proper 
value,  principal  value  and  eigenvalue  are  equivalent  terms.   The 
matrix  equation  that  must  be  solved  is  the  following: 

(B  -  AG_1)Q  =  0  (7) 

where  B  is  the  matrix  of  anisotropic  temperature  coefficients 

X  is  a  characteristic  value  of  the  equation 

G   is  the  reciprocal  space  metric  tensor  with  elements 

g..  =  b.  •  b.  where  the  b's  are  the  reciprocal  space  axis 
vectors  and 
Q  is  a  column  matrix  of  the  coefficients  of  the  characteristic 
vector.   This  is  a  vector  in  reciprocal  space. 

If  equation  (7)  is  premultiplied  by  G,  it  is  transformed  into  the 
standard  form  of  a  characteristic  value  problem. 

(GB  -  XI) Q  =  0  (8) 

where  I  is  the  identity  matrix  and 

G  is  the  real  space  metric  tensor  defined  in  equations  (3) 
and  (4). 

Equation  (8)  is  a  set  of  homogeneous  equations  of  order  3  which 
have  a  solution  if  the  determinant  of  the  coefficients  is  zero. 
Thus  the  characteristic  values  may  be  found  from  the  equation : 

49 


det  (GB  -  Al)  =  0  (9) 

Solution  of  this  equation  gives  three  values  of  A  which 
satisfy  equation  (8) .   In  order  to  have  a  meaningful  physical 
situation,  it  is  necessary  that  the  A's  which  are  the  roots  of 
equation  (9)  be  real  and  positive.   If  they  are  not,  the  set 
of  temperature  coefficients  is  said  to  be  non-positive  definite. 
This  is  sometimes  referred  to  as  complex  amplitudes  of  vibration 
or  negative  temperature  factors. 

Each  of  the  roots  of  equation  (9)  can  then  be  inserted  into  equat- 
ion (8)  and  the  resulting  equation  solved  for  the  characteristic 
vector  by  the  method  of  cof actors.  After  the  coefficients  of 
Q  are  determined,  they  may  be  used  to  form  the  reciprocal  space 
vector  q_  and  the  angle,  0^,  between  this  vector  q_  and  the  real 
space  axis  vector  a^  may  be  found  by  forming  the  dot  product 
of  the  two  vectors : 

%  •  £=  lajlsj  cos  0±         (10) 
This  equation  may  be  simplified  since 

%•£=%•    fajki  +  ^2-2  +  q3-3)  (11) 


and 


thus 


,  1  for  i  =  j  ,10v 

a.    .   b,   =  w.  .   =      _    ,.        ./.  (12) 

-d.       ^  ij  0  for  i  +  j 


a.  =  q,  (13) 


Equation  (12)  is  the  definition  of  the  reciprocal  space  axis 
vectors.  Equation  (10)  can  now  be  rewritten 

cos     Q±  =       qi  (14) 

IsJIal 

The  root-mean-square  amplitude  of  vibration  along  the  i'th  charact- 
eristic vector  is  called  y  and  is  found  from  the  i'th  root,  A  of 
equation  (9)  by  the  following  relationship: 

i±  =  CX± 

4.   Error  Analysis; 


y,  =   (A4/27T2)1/2  (15) 


If  y  is  a  computed  value  which  is  a  function  of  experimentally 

determined  variables  x, ,  x„,....x  such  that: 

12      n 

y=  y  (xj^,  x2 Xj^)  (16) 
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2 
then  the  variance  of  y,   denoted  by  a       is   found  from  the  equation 

n         n 

a       =     E          I  3y   9y  e              s                                      ,,_,. 

y          •    t      •    i  a     T~~  cov(x       x.)                                       (17) 

1=1      1=1  dx.dx.  11 

J              i     j  J 

where  cov(x  ,x.)  is  the  ij'th  element  in  the  variance-covariance 
matrix  of  the  x's. 

In  crystallographic  calculations,  the  variables  for  which  the 
variance-covariance  matrix  elements  are  known  include  the  unit 
cell  parameters,  the  atomic  positions  and  the  temperature  factor 
coefficients.   If  the  dependence  of  the  bond  distance  on  temper- 
ature factors  is  neglected,  as  is  done  in  this  program,  the  error 
in  the  bond  distance  will  depend  only  on  the  cell  parameters  and 
the  atomic  positions.   Similarly  the  ellipsoids  will  depend  only 
on  the  unit  cell  parameters  and  the  temperature  factor  coefficients. 
Thus,  for  crystallographic  problems,  equation  (17)  may  be  rewritten 
in  the  form: 


a2 


=  £    I   r-* %—      cov(a.,a.) 


y   j  i  ■  i   3aJ  9a. 
1=1  j=l   i    j 


(18) 


where  the  a_' s  are  the  unit  cell  parameters  and  the  z_'s  are  the 
atomic  parameters.   The  variance-covariance  matrix  is  the  inverse 
of  the  normal  equations  least  squares  matrix.   The  error  formula 
in  equation  (18)  is  a  "full  matrix"  error  formula.   The  so-called 
"diagonal"  form  includes  only  the  terms  in  equation  (18)  which 
have  i= j .   This  approximation  usually  gives  a  smaller  error.   All 
partial  derivatives  required  for  computing  the  errors  are  analyti- 
cally evaluated  in  this  program.   This  procedure  is  better  than 
numerical  differentiation  because  it  is  more  accurate,  easier  and 
faster  to  do  and  allows  singularities  in  the  derivatives  to  be 
detected  easier.   The  necessary  derivatives  are  evaluated  using 
the  chain  rule  of  derivatives. 

E.   OUTPUT 

1.   General  Output: 

The  first  output  consists  of  the  information  contained  on  the 
Title  Card.   Then  the  real  space  axis  lengths,  the  cosines  of  the 
real  space  angles ,  the  reciprocal  space  axis  lengths  and  the  cos- 
ines of  the  real  space  angles,  the  reciprocal  space  axis  lengths 
and  the  cosines  of  the  reciprocal  angles  are  printed.   The  real 
space  cell  volume  is  also  printed.   Then  the  lattice  type  and 
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centric  indicator  is  printed  followed  by  the  symmetry  operations. 
Next,  if  errors  are  to  be  calculated,  the  standard  deviations 
of  the  cell  parameters  and  the  cell  parameter  correlation 
matrix  are  printed.   Then  the  atomic  positions  and  temperature 
factors  are  printed.   If  ellipsoids  and  errors  are  to  be  cal- 
culated, the  standard  deviations  of  the  thermal  parameters  and 
temperature  factors  are  printed.   If  ellipsoids  and  errors  are 
to  be  calculated,  the  standard  deviations  of  the  atomic  positions 
and  the  correlations  matrix  for  the  positions  are  then  printed. 
This  completes  the  general  output. 

2.  Bond  Distance  and  Angle  Output: 

For  each  atom  in  the  asymmetric  unit ,  all  atoms  within  the  bond 
distance  limits  are  generated.   These  atoms  are  added  to  a  table 
of  generated  positions.   For  each  bond  distance  within  the  al- 
lowable range,  the  atom  identifier,  the  fractional  coordinates, 
the  bond  distance,  the  error,  if  calculated,  and  the  number  of 
this  atom  in  the  table  are  printed.   In  the  bond  angle  section 
of  the  program,  all  atoms  except  those  in  the  asymmetric  unit 
will  be  referenced  by  this  atom  number  rather  than  by  coordinates. 
This  number  makes  it  easy  to  determine  quickly  what  angle  has 
been  calculated. 

The  bond  angle  section  output  consists  of  the  atom  identifier 
and  coordinates  of  the  central  atom  which  is  one  of  those  in 
the  asymmetric  unit.   The  peripheral  atoms  are  each  referenced 
by  an  atom  identifier  and  an  atom  table  number.   The  bond  angle 
and  its  error,  if  calculated,  complete  the  output  for  each  angle. 

3.  Ellipsoid  of  Vibration  Output: 

For  each  anisotropic  atom,  the  equivalent  isotropic  B  and  its 
error  are  calculated  and  printed.   Then  for  each  of  the  char- 
acteristic values,  the  r.m.s.  amplitude  of  vibration  and  the 
angles  that  the  associated  characteristic  vector  makes  with  the 
real  cell  axes  as  well  as  the  errors  in  these  quantities  are 
printed.   No  output  occurs  for  isotropic  atoms. 

A.   Page  Estimates: 

The  main  program  will  produce  approximately  one  page  of  output 
plus  1/2  page  per  atom  if  errors  are  to  be  calculated.   The 
bond  distance  and  angle  subroutine  will  require  about  one  page 
for  each  60  distances  or  angles  calculated.   However,  the  user 
should  keep  the  bond  distance  limits  as  small  as  possible  since 
if  n  distances  are  within  the  limits,  the  number  of  angles  com- 
puted will  be  approximately  1/2  n  .   The  ellipsoids  of  vibra- 
tion subroutine  will  require  one  page  for  each  three  atoms. 
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PROGRAM  ERRORS  AND  ERROR  STOPS 

1.   There  are  two  kinds  of  errors  which  return  control  to  the 

monitor  through  a  STOP  statement.   If  a  symmetry  card  cannot 
be  interpreted  because  the  information  is  not  of  the  correct 
form  or  because  it  is  punched  in  the  wrong  columns,  the  pro- 
grams executes  a  STOP  11  statement  and  returns  to  the  monitor. 
The  card  in  error  will  be  the  last  symmetry  qard  listed  on  the 
output.   If  the  user  attempts  to  read  more  than  24  symmetry 
cards,  the  program  will  return  to  monitor  with  a  STOP  22  instruc- 
tion.  A  different  kind  of  error  can  occur  in  the  bond  distance 
calculation.   If  the  distance  limits  are  such  that  more  than 
1000  atom  positions  are  generated,  the  program  prints  a  dia- 
gnostic message  and  returns  to  the  main  program.   This  error 
does  not  terminate  the  program  and  ellipsoids  of  vibration  may 
be  calculated  correctly  after  this  condition  occurs.   However, 
bond  angles  will  not  be  calculated. 

COMMON  VARIABLES  AND  DESCRIPTION 

XYZ(3,60)  -  Fractional   coordinates  of  atoms  in  asymmetric  unit 

BETA(6,60)  -  Temperature  factor  coefficients 

A (3)  -  Real  cell  axis  lengths 

ANG(3)  -  Cosines  of  the  real  cell  angles 

AS (3)  -  Reciprocal  axis  lengths 

CS(3)  -  Cosines  of  the  reciprocal  cell  angles 

ATOM(2,60)  -  Atom  identifiers 

SIGA(6,6)  -  Variance-covariance  matrix  for  real  cell  parameters 

SIGB(6,6,60)  -  Variance-covariance  matrix  for  thermal  coefficients 

RS(3,3,48)  -  Rotational  matrix  part  of  symmetry  operator 

TS(3,48)  -  Translational  part  of  symmetry  operator 

TP(4,3)  -  Location  of  lattice  points  in  fractional  coordinates 

ITF(60)  -  Temperature  factor  type  indicator 

CORR(16290)  -  Variance-covariance  matrix  for  atomic  positions. 
This  is  stored  in  packed  form. 

XYZP(3,1000)  -  Fractional  coordinates  of  atoms  generated  in  the 
bond  distance  subprogram. 

IATP(IOOO)  -  Atom  identifier  of  atoms  generated  in  bond  dis- 
tance calculation. 

NA  -  Number  of  atoms  in  asymmetric  unit 

NSS  -  Number  of  symmetry  operations 

VOL  -  Volume  of  real  space  unit  cell 

DLIMIT(3)  -  Array  of  bond  distance  limits 

PI  -  The  value  of  the  constant  it 

G(3,3)  -  Real  space  metric  tensor 

NLP  -  Number  of  lattice  points  in  the  unit  cell 

IERR  -  Error  calculation  indicator 

TITLE (18)  -  Alphanumeric  title  information 
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H.   PROGRAM  NAMES  AND  DESCRIPTION 

1.  U  M  BADTEA: 

This  is  the  main  program  which  handles  the  input  and  calls 
the  bond  distance  and  ellipsoid  of  vibration  programs. 

2.  BODAN(IBOD): 

This  subprogram  calculates  bond  distances  and  errors  and  if 
the  argument  IBOD  is  equal  to  2 ,  calculates  the  bond  angles 
and  errors . 

3.  ELVIB: 

This  subprogram  calculates  the  ellipsoid  of  vibration  para- 
meters and  associated  errors. 

4.  VARIAN(A,B,C,D,E,II): 

This  subprogram  calculates  the  variance  of  a  computed  value 
where  the  partial  derivatives  of  the  value  with  respect  to 
the  cell  parameters  are  stored  in  array  B,  the  partial  deri- 
vatives with  respect  to  the  structure  parameters  are  stored 
in  array  C,  the  cell  convariance  matrix  is  stored  in  array 
D,  the  structure  parameters  covariance  matrix  is  in  array  E, 
the  argument  II  contains  the  atom  number  when  the  variance 
of  ellipsoid  parameters  is  calculated  otherwise  it  contains 
1  and  the  computed  variance  is  output  in  argument  A. 

5.  CELLD(B,A,AN,R,CS,VOL) : 

This  subprogram  computes  the  partial  derivatives  of  the  re- 
ciprocal cell  parameters  with  respect  to  the  real  cell 
parameters.   These  derivatives  are  used  in  computing  the 
ellipsoid  errors.   Array  B  contains  the  partial  derivatives 
after  execution.   Arrays  A,  AN,  R  and  CS  are  the  input  values 
of  the  real  cell  axis  lengths,  the  real  cell  cosines,  the 
reciprocal  axis  lengths  and  cosines  of  the  reciprocal  cell 
angles  respectively.   The  argument  VOL  is  the  input  value  of 
the  real  cell  volume. 

6.  DERDET(A,B,C)  : 

This  subprogram  calculates  the  partial  derivatives  of  matrix 
B  with  respect  to  the  thermal  coefficients  and  stores  them 
in  array  A.   Array  C  contains  an  intermediate  set  of  partial 
derivatives. 
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7.   C0RR1(SIGZ,NN,II,JJ,C0RR,KK,LL,MM) : 

This  subprogram  prepares  the  atomic  positional  correlation 
matrix  for  atoms  II  and  JJ  from  the  packed  matrix  CORR  of  rank 
MM  and  stores  it  at  the  KK,LL  position  of  the  array  SIGZ  which 
is  of  rank  NN. 

SYSTEM  REQUIREMENTS 

U  M  BADTEA  requires  two  tape  units,  a  standard  input  tape  and 
a  standard  output  tape.   All  input  is  BCD  from  Logical  Unit 
Number  5  while  output  is  BCD  to  Logical  Unit  Number  6.   The 
program  requires  approximately  32K  storage  locations.   If  the 
modified  ORFFE  tape  is  to  be  used,  it  is  read  from  logical 
unit  4. 
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Appendix  C 

LISTINGS  OF  PROGRAMS 

Fortran  listings  of  all  programs  and  subprograms  are  included 
herein.   Example  runstreams,  with  the  output  they  produce,  are  available 
on  request  from  either  author. 
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